Insights into the Maternal Ancestry of Côte d’Ivoire Honeybees Using the Intergenic Region COI-COII

Honeybee populations in Côte d’Ivoire have been previously identified as belonging to one subspecies, Apis mellifera scutellata, but other studies have since reported a mixed population consisting of A. m. adansonii and A. m. jemenitica. The population structure and the geographic distribution of honeybees in Côte d’Ivoire remain unclear. This study aimed to profile the population structure of honeybees and their biogeography in Côte d’Ivoire. A total of 33 honeybee colonies were sampled from 15 localities to investigate the maternal ancestry of indigenous honeybee populations using the DraI COI-COII mtDNA test. The results revealed that the honeybee population in Côte d’Ivoire is composed of African haplotypes, all belonging to the AI sublineage. Haplotypes A1 and A4 were recorded with five new sequence variants, including three types of haplotype A1 and two types of haplotype A4. The A1e variant was the most frequent in the A. m. adansonii distributional area. The distribution of the haplotype variants was correlated with the climate pattern in Côte d’Ivoire. This is the first study in Côte d’Ivoire that gives insights into the biogeography and mitotype structure of the local honeybee populations.


Introduction
Based on morphometrics and multivariate analyses [1], the honeybee species Apis mellifera L. has been split into four lineages. Lineage A is spread from the north to the south of Africa, lineages M and C are distributed in Europe and lineage O in the Middle East [1][2][3]. Recently, the DraI mitochondrial DNA test was used in the identification of maternal honeybee populations [3][4][5][6][7][8][9]. Using this method, four mitochondrial DNA (mtDNA) lineages have been identified, including lineages A, M, C, and Y [2,7,10]. Honeybees from the M and C lineages have been more extensively studied for their taxonomy, biodiversity, and geographic distribution than their counterparts, the A and Y lineages. The genetic diversity of honeybees in the African continent is still understudied [3,11,12], making genetic conservation programs difficult in Africa [12]. To gain further insights into the genetic diversity of indigenous honeybee populations, more studies are needed to improve the availability of reference data in different African regions. This information could contribute towards monitoring the endemic populations, improving honeybee services, and designing strategies for ecological conservation of the local population. Several studies have been done to discriminate the honeybee population in Africa [1,3,[13][14][15][16][17][18][19][20][21][22]. Using morphometric traits, eleven A. mellifera L. subspecies have been taxonomically recognized in the African continent, i.e., A. m. lamarckii [14], A. m. jemenitica [15], A. m. intermissa [16], A. m. sahariensis [17], A. m. unicolor [18], A. m. scutellata [19], A. m. capensis [20], A. m. adansonii [18], A. m. monticola, A. m. litorea [21], and A. m. sinensis [22].
In West Africa, particularly in Côte d'Ivoire, the biodiversity of honeybee populations has been studied with sample sizes that are low for the size of the region (322,462 km 2 ). Classification of local honeybee populations in Côte d'Ivoire has, therefore, been controversial. Based on the morphometric characterization, the honeybee populations in Côte d'Ivoire were initially claimed to be A. m. scutellata [19]. Later, Ruttner [1] recognized these populations as belonging to A. m. adansonii in the tropical dry and subequatorial climates. Similarly, according to the morphometric study carried out by Radloff et al. [13], honeybee populations from Côte d'Ivoire were identified as A. m. adansonii in the tropical dry and tropical humid climates, with some A. m. jemenitica and A. m. adansonii hybrid populations.
Mitochondrial DNA variation in the cytochrome oxidase subunit I-II (COI-COII) intergenic region has been used to distinguish lineages and to refine the classification of A. mellifera L. [2,[23][24][25][26]. Mitochondrial DNA markers were employed to demonstrate that the honeybee populations from Côte d'Ivoire belong to lineage A [3,27] [28], and Z (haplotype Z 1 -Z 4 , Z 7 ) [29,30]. These sublineages are usually differentiated by the presence or absence of an additional DraI site (TTTAAA) and a deletion at the 3' end of the P element [2,3]. Two forms of the P element, P 0 and P 1 , are typical of lineage A. The P 1 form is characterized by a 15-bp deletion at the 3' end of the P element, whereas P 0 does not exhibit any large deletion. The P 0 form is carried by sublineages A I , A II , and Z [3,29,30], whereas P 1 is carried by sublineage A III [28]. Sublineage A II is differentiated from sublineage A I by the absence of the DraI site at the 5' end of the first Q element, whereas sublineage Z has an additional DraI site in the middle of the first Q element. The classification into sublineages and the haplotype geographical distribution patterns of the honeybee populations of Côte d'Ivoire are yet to be clarified. Getting an accurate view of the genetic diversity could allow targeting of the sensitive honeybee populations and conservation of their biodiversity.
In this study, we provide for the first time the matriline structure and biogeographic distribution of the honeybee populations in Côte d'Ivoire using the highly polymorphic intergenic COI-COII region of the mtDNA.

Study Area
The study was conducted in Côte d'Ivoire, a West African country geographically located at 5 • 18 34" N and 4 • 00 45" W, between the Sahara Desert and the Atlantic Ocean. The country has four climate types from the north to the south (Figure 1), i.e., tropical dry climate, tropical humid climate, subequatorial climate, and a mountain climate that is found on the western side (Table 1).  Table 1.

Collection of Honeybees
Between May and July, 2015, sixty adult worker honeybees were collected from each of 33 colonies from 15 localities in Côte d'Ivoire (Table 1). Ten colonies were sampled from the tropical dry, tropical humid, and subequatorial climates each, and three colonies were sampled from the mountain climate (Table 1). Samples were taken from feral swarms collected by the traditional beekeepers in the countryside. All samples were preserved in 90% ethanol and were transferred under cool conditions to the Institute of Apicultural Research (IAR), Chinese Academy of Agricultural Sciences, Beijing, China, for their molecular characterization. The low sample size collected in Côte d'Ivoire was mainly due to the quasi absence of beekeeping in several regions. Moreover, in the selected regions, beekeepers (only 18) held less than five hives. Each sampling locality is presented with a black dot. The 15 sampling localities were spread over the 4 climate types (tropical dry and humid climates, subequatorial dry, and mountain climates) and the 2 biome types (Savannah and Forest). The numbers in the map represent a site with correspondence names and coordinates shown in Table 1.

Collection of Honeybees
Between May and July, 2015, sixty adult worker honeybees were collected from each of 33 colonies from 15 localities in Côte d'Ivoire (Table 1). Ten colonies were sampled from the tropical dry, tropical humid, and subequatorial climates each, and three colonies were sampled from the mountain climate (Table 1). Samples were taken from feral swarms collected by the traditional beekeepers in the countryside. All samples were preserved in 90% ethanol and were transferred under cool conditions to the Institute of Apicultural Research (IAR), Chinese Academy of Agricultural Sciences, Beijing, China, for their molecular characterization. The low sample size collected in Côte d'Ivoire was mainly due to the quasi absence of beekeeping in several regions. Moreover, in the selected regions, beekeepers (only 18) held less than five hives.

DNA Extraction, Sequencing, and DraI Test
One worker honeybee from each sampled colony was subjected to DNA extraction. The DNA was extracted from the individual honeybee thoraces by using an E.Z.N.A ® Tissue DNA Kit (Omega Bio-Tek, Doraville, GA, USA) according to the manufacturer's instructions.
PCR amplification of the COI-COII intergenic region was carried out according to a protocol detailed previously [4]. The size of the amplified DNA amplicon was determined by running 10 µL of the PCR-amplified products for 20 min on 1.0% agarose gel using gel electrophoresis. After gel migration, it was visualized and photographed using a UV-equipped gel documentation system (Bio-Rad laboratories, 6000; Biorad, Hercules, CA, USA). About 20 µL of PCR products from each sample were sent to Sangon Biotech (Beijing, China) for purification and direct Sanger sequencing in both directions. To confirm the result quality of the sequenced DNA, we repeated the same analysis independently with two more individuals from each colony, yielding a total of three sequences per colony.
To conduct the restriction fragment length polymorphism (RFLP) analysis, we used DraI (Boehringer Manheim) to digest 500 ng of the PCR products from each colony individually. The digested products were then run on 2.0% Metaphor agarose gels prepared in 1 × TBE at 10 V/cm to carry out the electrophoretic analysis. The nucleotide bands were visualized under a UV transilluminator [31].

Data Analysis
Before data analysis, we used the CodonCode aligner (www.codoncode.com) to clean the sequences. Then, the mtDNA sequences were aligned using MEGA 5.04 software [32]. Comparisons of the sequences were conducted with the Basic Local Alignment Search Tool (BLAST ® ) by searching the most relevant DNA sequences available on the GenBank ® web portal (http://blast.ncbi.nlm.nih. gov/Blast.cgi). Newly described sequences were submitted to the GenBank NCBI database under the accession numbers MF984182, MF984186, MH152663, MH152664, and MH152665. Similarity among the COI-COII haplotypes was investigated using the PopART version 1.6 (Population Analysis with Reticulate Trees) software (http://popart.otago.ac.nz). The classical phylogenetic tree, including lineages M, C, and Y, and sub-lineages A I , A II , A III , Z, was built. We trimmed all nucleotide sequences to the same length. Then, the data was aligned on the online service of www.ebi.ac.uk. Using the software Jalview version 2.10.3 [33], the parsimony tree of mtDNA COI-COII was generated by average distance with the percentage identity method (PID).

Results
Honeybee populations from Côte d'Ivoire were composed of two mtDNA sequences, P 0 Q and P 0 QQ, corresponding to the sequence lengths of 545 and 737 bp, respectively. P 0 Q and P 0 QQ types were described by the haplotype A 1 [2] and A 4 respectively. While P 0 Q was the most recorded with an overall presence in 31 out of 33 sampled honeybee colonies, P 0 QQ was found in 2 out of 33. Among the 33 colonies examined, we recorded four variants of haplotype A 1 , i.e., A 1e [34], A 1s , A 1q , and A 1r , and two variants of haplotype A 4 , i.e., A 4o and A 4r . In this study, the haplotype variants A 1s , A 1q , A 1r , A 4o , and A 4r were newly described. The sequence A 1e was the most dominant (23 out of the 33) in the samples.
The sequences A 1 , A 1e , A 1s , A 1q , A 1r , A 4o , and A 4r comprised one haplogroup by network analysis (Figure 2). The sequence A 1e represents the internal component of the haplogroup. The sequences A 1 , A 1s , A 1q , and A 1r diverged from the dominant sequence A 1e by only one base ( Figure 3) and the sequences A 4o and A 4r were separated from A 1e by two bases.

of 11
aplotype network displaying the relationship between the mtDNA sequences of Apis mellifera L. samples collected from different localities in Cote represents the number of haplotype copies recorded in the dataset. Each hatch mark represents a nucleotide change.   Our results clearly show that Côte d'Ivoire honeybee colonies belong to the A lineage and the A I sublineage. The P 0 Q sequence pattern of A 1e , A 1s , A 1q , and A 1r and the P 0 QQ sequence pattern of A 4o and A 4r from the collected data were genetically close to the sequence pattern of the honeybee A. m. scutellata mtDNA of haplotype A 4 recorded in South Africa [5] and A. m. adansonii matrilines of haplotype A 1 recorded in Zambia [3] (Figure 4).
Insects 2019, 10 8 of 11 Our results clearly show that Côte d'Ivoire honeybee colonies belong to the A lineage and the AI sublineage. The P0Q sequence pattern of A1e, A1s, A1q, and A1r and the P0QQ sequence pattern of A4o and A4r from the collected data were genetically close to the sequence pattern of the honeybee A. m. scutellata mtDNA of haplotype A4 recorded in South Africa [5] and A. m. adansonii matrilines of haplotype A1 recorded in Zambia [3] (Figure 4). The sequence A1e was distributed across all climatic types (i.e., dry tropical, humid tropical, mountain, and subequatorial climates) in both savannah and forest biotypes ( Figure 5) with a frequency of six in the dry tropical climate, eight in the humid tropical climate, one in the mountain climate, and eight in the subequatorial climate. We found that A1r was only distributed in the dry tropical climate in the savannah, and A1q and A1s were recorded only in the humid tropical climate, in the savannah biome. Haplotype A1 was recorded in the subequatorial climate in the forest biotype, while the haplotype A4 and its variants A4o and A4r were predominant in the northern part of the country in the tropical dry climate of the savannah biotype. With the exception of A1e, there was a clear co-segregation of haplotype variants and climatic types. The sequence A 1e was distributed across all climatic types (i.e., dry tropical, humid tropical, mountain, and subequatorial climates) in both savannah and forest biotypes ( Figure 5) with a frequency of six in the dry tropical climate, eight in the humid tropical climate, one in the mountain climate, and eight in the subequatorial climate. We found that A 1r was only distributed in the dry tropical climate in the savannah, and A 1q and A 1s were recorded only in the humid tropical climate, in the savannah biome. Haplotype A 1 was recorded in the subequatorial climate in the forest biotype, while the haplotype A 4 and its variants A 4o and A 4r were predominant in the northern part of the country in the tropical dry climate of the savannah biotype. With the exception of A 1e , there was a clear co-segregation of haplotype variants and climatic types.

Discussion
Among the sequences of haplotype A1 and A4 (i.e., A1e, A1s, A1q, A1r, A4o, and A4r) recorded in this study from the Côte d'Ivoire, three variants of haplotype A1, namely A1s, A1q, and A1r, and two variants of haplotype A4, namely A4o and A4r, were not yet recorded and are described in this study. The haplotype A1 and its variant A1e were previously reported by Chávez-Galarza et al. [2] from Iberia and by Szalanski and Magnus [34] from Washington and Iron Counties, Utah, USA, in NCBI GenBank using mtDNA markers. All of them belonged to the evolutionary lineage A, particularly to the AI sublineage, indicating that honeybees from Côte d'Ivoire belong exclusively to the African lineage [3,13]. The hybridization with A. m. jemenitica from lineage Y, as reported earlier, may occur, between A. m. adansonii queen and A. m. jemenitica drone. Therefore, haplotypes of Y lineage ancestry (A. m. jemenitica) could not be found in our samples. The absence of Y lineage in our study does not mean that they are not present at all in Côte d'Ivoire. The rarity of A. m. jemenitica could be due to an inadequate environment, mainly due to low altitudes (ranging from 20 to 385 m) within the country. This result corroborates well with the previous work of El-Niweiri and Moritz [35] in Sudan, who reported that the populations of A. m. jemenitica are scarce at altitudes below 500 m.
The absence of honeybees from the M and C mtDNA lineages could be explained by the fact that these lineages are not endemic to West Africa [1,3]. In fact, the importation program of honeybees from the M lineage and Z sublineage initiated by the Government of Côte d'Ivoire in 1980 to improve the production of honey in the country had failed. The unsuccessful breeding of European honeybees in Côte d'Ivoire was due to the aggressiveness of the endemic honeybees against the imported honeybees and also to their poor adaptation to the new environment [36]. These factors significantly reduced the number of worker honeybees imported until the end of the honeybee importation program in Côte d'Ivoire in 1983. Therefore, the probability that these M and C lineages, and Z sublineage, were collected in our samples was very low.
The haplotype variants A4 occurred in the tropical dry climate in the ecological region of A. m. scutellata [35,36]. Haplotype A1 and its variants (A1e, A1s, A1q, and A1r) were recorded in the distributional area of A. m. adansonii [1,15,35]. The geographical distribution of the sequences recorded in this study is marked by a co-segregation of haplotypes and climate types. The genetic diversity of honeybee populations found in this study, therefore, might be a result of a long period of ecological adaptation to the environment. Thus, it is crucial to conserve the diversity of honeybee populations in Côte d'Ivoire.

Discussion
Among the sequences of haplotype A 1 and A 4 (i.e., A 1e , A 1s , A 1q , A 1r , A 4o , and A 4r ) recorded in this study from the Côte d'Ivoire, three variants of haplotype A 1 , namely A 1s , A 1q , and A 1r , and two variants of haplotype A 4 , namely A 4o and A 4r , were not yet recorded and are described in this study. The haplotype A 1 and its variant A 1e were previously reported by Chávez-Galarza et al. [2] from Iberia and by Szalanski and Magnus [34] from Washington and Iron Counties, Utah, USA, in NCBI GenBank using mtDNA markers. All of them belonged to the evolutionary lineage A, particularly to the A I sublineage, indicating that honeybees from Côte d'Ivoire belong exclusively to the African lineage [3,13]. The hybridization with A. m. jemenitica from lineage Y, as reported earlier, may occur, between A. m. adansonii queen and A. m. jemenitica drone. Therefore, haplotypes of Y lineage ancestry (A. m. jemenitica) could not be found in our samples. The absence of Y lineage in our study does not mean that they are not present at all in Côte d'Ivoire. The rarity of A. m. jemenitica could be due to an inadequate environment, mainly due to low altitudes (ranging from 20 to 385 m) within the country. This result corroborates well with the previous work of El-Niweiri and Moritz [35] in Sudan, who reported that the populations of A. m. jemenitica are scarce at altitudes below 500 m.
The absence of honeybees from the M and C mtDNA lineages could be explained by the fact that these lineages are not endemic to West Africa [1,3]. In fact, the importation program of honeybees from the M lineage and Z sublineage initiated by the Government of Côte d'Ivoire in 1980 to improve the production of honey in the country had failed. The unsuccessful breeding of European honeybees in Côte d'Ivoire was due to the aggressiveness of the endemic honeybees against the imported honeybees and also to their poor adaptation to the new environment [36]. These factors significantly reduced the number of worker honeybees imported until the end of the honeybee importation program in Côte d'Ivoire in 1983. Therefore, the probability that these M and C lineages, and Z sublineage, were collected in our samples was very low.
The haplotype variants A 4 occurred in the tropical dry climate in the ecological region of A. m. scutellata [35,36]. Haplotype A 1 and its variants (A 1e , A 1s , A 1q , and A 1r ) were recorded in the distributional area of A. m. adansonii [1,15,35]. The geographical distribution of the sequences recorded in this study is marked by a co-segregation of haplotypes and climate types. The genetic diversity of honeybee populations found in this study, therefore, might be a result of a long period of ecological adaptation to the environment. Thus, it is crucial to conserve the diversity of honeybee populations in Côte d'Ivoire.

Conclusions
According to our analyses of the intergenic COI-COII region of mtDNA, honeybee populations in Côte d'Ivoire exhibited the occurrence of two haplotypes, A 1 and A 4 , with different variants i.e., A 1e , A 1s , A 1q , A 1r , A 4o , and A 4r . We did not find any local evidence of the introgression of imported honeybees from M and C lineages and Z sublineage. The haplotypes A 1 and A 4 occurred in the ecological regions of A. m. adansonii and A. m. scutellata, respectively. However, we were unable to fully infer the taxonomy of our honeybee colonies. Therefore, we recommend that future research should combine both morphometric and nuclear DNA analyses such as microsatellite markers for the accurate identification of the local honeybees in Côte d'Ivoire. Nevertheless, a larger sampling area is needed to be prospected to get a deeper insight into the biogeography of A. mellifera L. in Côte d'Ivoire.