Origin, evolution and conservation of the honey bees from La Palma Island (Canary Islands): molecular and morphological data

Genetic studies have shown a significant hybridization of local populations of Apis mellifera in the Canary Islands, as a consequence of queen importation. Since La Palma Island was a probable exception to this hybridization rule, a local honey bee Conservation Project was developed. As a first step, the genetic characterization of the population was performed with the scope of generating fundamental knowledge to correctly manage honey bee genetic resources. In this sense, the origin of the Canarian honey bee is at issue, since diverse morphological and genetic studies have shown it to be either African or European. In our study, 499 colonies from La Palma were analyzed using mitochondrial DNA, microsatellites and wing geometric morphometric data. Hybridization from the C evolutionary lineage was observed mainly at mitochondrial level, but this hybridization was focused on a restricted area and hybrid colonies were detected for replacement. All genetic and geometric morphometric results showed the Canarian honey bee clustered with Western European M branch, specifically with A. m. iberiensis subspecies, far from African ones. We concluded that the present Canarian honey bee gene pool would have originated most likely from human introductions from Portugal dating back to the conquest of the islands in the sixteenth century, although an ancestral natural colonization from Iberia cannot be ruled out. Results also showed that Canarian honey bees are genetically differentiated from A. m. iberiensis, probably due to micro-evolutionary processes, such as founder effects, bottleneck and local adaptation. Therefore, we propose the classification of these local populations as a differentiated ecotype of A. m. iberiensis. This study highlights the idiosyncrasy of the Canarian honey bee and the need to protect it, as it is a population that extends considerably the genetic diversity of M branch.

Based on a "local solutions to global problems" proposition, loss of biodiversity across the planet has engendered a special interest in recovery and conservation of local populations of animals, and A. mellifera is no exception. Given the regression of local honey bee populations in some areas, several local initiatives have emerged, including projects for the conservation of local honey bees in several European countries. The Canary Islands are one such example. De La Rúa et al. (1998) reported high levels of hybridization of the Canarian honey bee with foreign honey bees from the C evolutionary branch. Interestingly, hybridization was not homogeneously distributed among the Canary Islands while it was as high as 39% in the island of Tenerife, the island of La Palma showed no signs of hybridization (De La Rúa, Serrano, Galian, & Moritz, 2001;De La Rúa et al., 1998). This finding led to the creation of a local Honey Bee Conservation Project on La Palma, the most western island of the Canarian archipelago. With an area of 706 km 2 , it is 450 km from the African coast and more than 70 km from any other Canarian island. In 2002, La Palma has been designated as a Biosphere Reserve. Current law protects local populations and honey bee importation to La Palma is forbidden (Government of Canary Islands, 2001: BOC 49 of 20/4/ 2001. In addition, the conservation strategy includes an Area of Special Protection in the northeast of the island, managed by the local Honey Bee Conservation Project, establishing a natural mating area favored by a natural barrier that allows saturation with local drones. The genetic characterization of populations is a fundamental element for the correct management of animal genetic resources in a viable conservation plan, but the evolutionary history and taxonomic classification of the Canarian honey bee are not clearly resolved. Morphological studies have shown similarity between the bee populations from the Canary Islands and those from the Iberian Peninsula (Padilla, Hernández, & Reyes, 2001;Padilla et al., 1998;Ruttner, 1975), a result consistent with historical and cultural evidence. On the contrary, studies with genetic markers have reported the singularity of the Canarian honey bee populations and, consistent with their geographical proximity, these studies show a closer relationship with Moroccan bees (A branch) than with their Iberian counterparts (M branch) (De La Rúa et al., 1998, 2001. Excluding hybridization processes with C lineage, only mitochondrial haplotypes from A (African) branch have been detected, and some novel and distinct haplotypes were described by De la Rúa et al. (1998): haplotypes A15 and A14 (included later in A III sublineage by Franck et al., 2001). The more frequent Iberian A branch haplotypes were absent in the Canarian honey bee. All these results, along with the genetic similarity between Canarian honey bees and those from Madeira and Azores (De La Rúa, Galian, Pedersen, & Serrano, 2006) led De la Rúa et al. (1998Rúa et al. ( , 2001Rúa et al. ( , 2006 to suggest that the honey bees from Macaronesia were originally derived from African populations, and that thereafter they underwent genetic differentiation due to isolation from continental populations. Recent studies (Pinto, Muñoz, Chávez-Galarza, & De La Rúa, 2012;Pinto et al., 2013) have reported a high amount of novel A III sublineage haplotypes in Portugal, which opens the possibility of an Iberian origin of the Canarian A III sublineage haplotypes.
The goal of the present study was twofold: first, to increase the knowledge about the evolutionary history of the honey bees from the Canary Islands, analyzing both morphological and molecular data, and second, to provide information and practical genetic framework for animal management in the conservation program already underway on La Palma. To address these genetic issues, we screened for hybridization and analyzed genetic diversity and differentiation on a large sample of honey bee populations from La Palma. For comparatives purposes, we included European and African samples and performed phylogenetic analyses of populations based on morphological and molecular data.

DNA extraction and analysis
Mitochondrial DNA DNA was extracted from homogenized head and thorax of a single worker per colony using either a phenolchloroform method or a chelex-based protocol (Walsh, Metzqer, & Higuchi, 1991). Four hundred and eightyeight colonies from La Palma were analyzed for the COI-COII intergenic region (Supplementary Table 1), as described in Garnery, Solignac, Celebrano, and Cornuet (1993). This region is located between the cytochrome oxidase I (COI) and the cytochrome oxidase II (COII), and comprises the tRNA leu gene, P or P 0 and Q(s) sequences, and the 5´end of the COII gene (Garnery, Cornuet, & Solignac, 1992;Garnery, Mosshine, Oldroyd, & Cornuet, 1995). After amplification, fragment length was determined by electrophoresis in 1.4% agarose gel. The DraI restriction pattern was identified by 10% PAGE stained with ethidium bromide, and mitochondrial lineage and haplotype were determined. The COI-COII intergenic region presents a mix of single base mutations, large deletions, insertions, and/or duplications that provide limited insight into phylogenetic relationships of these populations. Accordingly, to study phylogenetic relationships of populations, a 715 bp fragment of the cytochrome b (Cytb) mitochondrial gene was sequenced in a total of 43 colonies: 15 from La Palma (all of them with COI-COII A lineage haplotypes) ( & Crozier, 1993). Amplification was carried out using the CB1 primer (5´-TATGTACTAC-CATGAGGACAAATATC-3´; Jermiin & Crozier, 1994), and TSER primer (5´-ACTTATTCAAGTTCATTAACT-3´) with the COI-COII intergenic region amplification protocol (Garnery et al., 1993) but adjusting annealing T˚to 49˚C. Amplification products were purified by the purification kit QIAquick (QIAGEN). The sequencing was carried out by ABI PRISM™ BIGDYE v3.1 â Terminator Sequencing Reaction â (Applied Biosystems) and an automatic sequencer ABI PRISM 3100 Genetic Analyzer. The sequences were edited using BioEdit v7.0.9 software (Hall, 1999). Alignment was achieved using ClustalW2 software (Larkin et al., 2007).

Geometric morphometrics
Geometric morphometrics uses landmark coordinates to achieve a better description of biological structures with the benefits of a higher statistical power than traditional morphometrics (Adams, Rohlf, & Slice, 2004;Bookstein, 1991). Applied to the honey bee wing shape, geometric morphometrics has proven useful for the analysis of wing asymmetry (Smith, Crespi, & Bookstein, Table 1. Haplotype diversity for COI-COII, expected heterozygosity for 10 microsatellites loci, and introgression with C lineage for mitochondrial and nuclear data (weighted using registered colonies), in sampled populations (Figure 1) 1997), the differentiation of A. mellifera subspecies (Francoy et al., 2008(Francoy et al., , 2009Kandemir, Ozkan, & Fuchs, 2011;Miguel et al., 2011;Tofilski, 2008), and the semiautomatic identification of bee workers using wing shape landmark coordinates (Baylac et al., 2008). In the present study, we applied geometric morphometrics to the analysis of wing shape described by the coordinates of 19 landmarks of the forewing already selected by Smith et al. (1997) and applied by Miguel et al. (2011). The used sequence and numbering of points is the recommended by Meixner et al. (2013). Temporary preparations of the right wings were done in distilled water. After careful removal of the excess water, digitized images were taken using a video camera connected to a framegrabber. Measurement TV software (Updegraff, 1990) was used to measure landmark coordinates. Each wing was measured twice and averages were used to minimize measurement errors. Wing terminology followed Michener (2000).

Statistical analyses
Mitochondrial DNA Haplotype diversity and standard error for mitochondrial DNA were calculated according to Nei and Tajima (1981). The phylogenetic trees of the 44 Cytb sequences were constructed using the PAUP 4.0.1.0 software (Swofford, 1996). The neighbor-joining algorithm (NJ) was carried out with the HKY85 substitution model (Hasegawa, Kishino, & Yano, 1985) and bootstrap values were calculated from 1000 replications. Maximum parsimony analysis was performed with random sequence addition (100 replicate heuristic searches), after which a consensus tree was done. Prior to maximum likelihood (ML) analyses, the MODELTEST 3.7 program (Posada & Crandall, 1998) with the Akaike information criterion, was used to choose the DNA substitution model best fitting our data. ML searches were conducted using the chosen model. Parsimony and likelihood bootstrap confidence values were calculated from 1000 to 100 replications, respectively, and tree bisection and reconnection (TBR) branch swapping heuristic search algorithm were employed. The Tajima neutrality test and divergence estimation were performed using DnaSp 3.99 software (Rozas, Sanchez-DelBarrio, Messeguer, & Rozas, 2003). Additionally, a median-joining network of 44 Cytb sequences was constructed using the software Network 3.1.1.1 (Bandelt, Forster, & Rö hl, 1999).
major, A. m. sahariensis, and A. m. intermissa from Morocco (Franck, Garnery, Solignac, & Cornuet, 1998), A. m. ligustica from Italy, A. m. macedonica from Greece (Estoup, Garnery, Solignac, & Cornuet, 1995), A. m. carnica from Germany and A. m. caucasica from Georgia (L. Garnery unpublished data). Heterozygosity was calculated according to Nei (1978). To estimate the nuclear hybridization of La Palma honey bees with C branch, Introgression Rate (IR) formula was applied using the frequencies of diagnostic alleles (Garnery et al., 1998b) for 8 of the 10 microsatellite loci set, because diagnostic alleles were not described for B124 and Ap55 microsatellites. The exact test for Hardy-Weinberg equilibrium was performed for populations of La Palma, using the GENEPOP package version 4.0 (Rousset, 2008). F ST and F IS were calculated according to Weir and Cockerham (1984). The confidence intervals were determined by jackknife and statistical significance applying 15,000 permutations, all obtained by FSTAT software (Goudet, 2001). When necessary, Bonferroni's correction was applied (Weir, 1996): probability of p < .006 and p < .001 was considered significant for 10 and 6 microsatellite loci, respectively. Population structure was studied using the STRUCTURE software (version 2.3.3; Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000). In this analysis, a Bayesian-based approach was used to search for the occurrence of independent genetic groups (K) with individuals allowed being products of admixture between two or more of the populations. The method employs a Markov chain Monte Carlo algorithm to estimate: (i) the allele frequencies for each K, and (ii) the proportion of the genome of each individual derived from each population (genetic component). To ensure that an uneven sample did not affect the results, and since there were only around 30 individuals of the subspecies of the A, C, and O evolutionary branches (Estoup et al., 1995;Franck et al., 1998), only one population of each of the subspecies of the M branch (selected at random) was included in the analysis. Thirty individuals from La Palma were likewise selected at random. A total of 20,000 Markov chain Monte Carlo iterations, after a burn-in of 20,000 iterations, were run 10 times for each number of genetic clusters (K). We chose the option of correlated allele frequencies between populations as recommended by Falush et al. (2003). The method proposed by Evanno, Regnaut, and Goudet (2005) was used to identify the best value of K. This method calculates K, an ad hoc statistic that is based on the rate of change in the log probability of data between successive K values. We used CLUMPP v1.1.2 software (Jakobsson & Rosenberg, 2007) to determine the optimal assignation of clusters for the analyzed individuals, maximizing similarity between the different STRUCTURE replications. The resulting clusters were visualized using DISTRUCT v1.1 software (Rosenberg, 2004). D A genetic distance for microsatellites (Nei, Tajima, & Tateno, 1983), NJ trees, and bootstrap procedures (1000 iterations) were performed using Populations 1.2.28 software (Langella, 1999).

Geometric morphometric data
Raw landmark coordinates were superimposed using a generalized least-squares procrustes superimposition algorithm (Rohlf & Slice, 1990), and the superimposed coordinates were projected onto the tangent shape space (Rohlf, 1999), thus providing the shape parameters used in the multivariate analyses. Shape differences between localities were investigated using canonical variate analyses (CVA). Wing deformations along factorial axes were estimated by the multivariate regression (Monteiro, 1999) and two extreme shapes were then calculated, one for each axis extremity. All correct classification rates were estimated using leave-one-out cross-validations (looCV), which provide lower but unbiased estimates (Ripley, 1996). Finally, a D2 Mahalanobis distance matrix (Mahalanobis, 1936) between populations was built. Procrustes superimposition and statistical analyses were done using geometric morphometric functions for MATLAB devised by one of us (M.B.).

Mitochondrial DNA
The COI-COII mtDNA haplotype distribution among honey bee populations from La Palma is shown in Figure 1 and in Supplementary Table 2. Haplotypes belonging to the A (African) evolutionary lineage were largely predominant. A15 was the most common haplotype (228 of 488, 47%). A15 is included in A III sublineage (Franck et al., 2001), which presents a 17 bp deletion in P 1 sequence as described by De la Rúa et al. (1998).
Haplotypes A11 (24 of 488, 4.9%), A14 (27 of 488, 5.5%), and A20 (1 of 488), which also belong to A III sublineage, were likewise detected in samples of La Palma. Another common haplotype found in La Palma was A1 (36%), classified in A I sublineage (Franck et al., 2001). COI-COII haplotype diversity values of La Palma (Table 1), with a mean value of .6456, were significantly lower than Iberian (.8107; Miguel et al., 2007) and African (.7164;Franck et al., 1998) ones (p < 0,001; Student t-test). European M lineage haplotypes were absent. On the contrary, the weighted mean frequency of the C lineage haplotypes in La Palma was 10.7% (Table 1). However, this haplotype was restricted to the central/western part (Figure 1), with frequencies reaching as high as 54.5% in El Paso. The C lineage haplotypes were not detected in the Area of Special Protection. Analysis of the Cytb gene (GenBank accession numbers of sequences are given in Supplementary Table 3) in 44 individuals identified 15 different haplotypes (H1-H15). D value from Tajima neutrality test was not significantly different from zero (D = −1.274; p > .1), i.e., sequence variation adjusted to a neutral evolutionary model. Maximum likelihood analysis of Cytb sequences (Figure 2) showed that all the mitochondrial sequences of A lineage clustered together, apart from M and C lineage sequences. Within the A lineage, La Palma and Iberian haplotypes grouped together and separated from all but one Moroccan ones. When COI-COII haplotypes were translated to the phylogeny of Cytb sequences, no direct relationship was observed between the two markers; six out of the 12 analyzed COI-COII haplotypes of A lineage (A1, A2, A3, A8, A14 and A15) were associated with more than one Cytb haplotype, while H07 Cytb haplotype included sequences of 11 COI-COII haplotypes (Figure 2). In the median-joining network built using the Cytb sequences (Figure 3), a large phylogenetic distance among A1 haplotypes from Africa (H12 and H13) and those from Iberia and La Palma (H07 and H08) was shown.

Microsatellite loci
Among the 12 populations of La Palma, 10 were found to be in Hardy-Weinberg equilibrium after applying Bonferroni correction. The two populations that were not in Hardy-Weinberg equilibrium (Barlovento and San Andrés) showed an excess of homozygotes, with positive and statistically significant F IS values (Barlovento, F IS = .061 ± .060; San Andrés, F IS = .028 ± .039). The number of detected alleles per microsatellite ranged from 3 (A28) to 10 (B124) (Supplementary Table 4). The highest heterozygosity values (above .30) were found for A113, Ap66, Ap81, and B124 loci, and the lowest (below .025) for A24, A28, and A88. Average heterozygosity values of the populations varied from H e = .175 in San Andrés to H e = .226 in Puntagorda (Table 1). Mean heterozygosity of 6 microsatellites (A7, A28, A113, B124, A24, A88) reflected that the diversity of La Palma (H e = .1495) was significantly lower than that observed in the Iberian Peninsula (H e = .3659; p < .001; Student t-test) and in northern Africa (H e = .7539; p < .001; Student t-test).
Populations were analyzed using STRUCTURE software to study the genetic relationship between La Palma honey bees with A. mellifera subspecies, and to identify the genetic contribution of European and African subspecies to La Palma honey bees. The best adjustment of the data after applying the procedure proposed by Evanno et al. (2005) was 3 clusters (Figure 4(A) and Supplementary Figure 1(A)). La Palma honey bee was grouped with M branch subspecies. The three African subspecies were grouped in a second cluster, and the subspecies of the C and O evolutionary branches appeared grouped in a third cluster. The STRUCTURE software also allows the individuals to be classified into the different clusters obtained. Thus, to balance the sample size in the analysis, the 489 La Palma individuals were divided into 12 groups of approximately 40 individuals, and these 12 groups were analyzed together with the rest of subspecies (results not shown). All La Palma individuals were classified in the M cluster, together with A. m. iberiensis and A. m. mellifera. Only 17 individuals showed admixture levels over 5%. The average admixture percentage was .8% with the African cluster and .6% with the C+O cluster. It should be stressed that 7.5% of the genome of the Lisbon population (A. m. iberiensis) presented a common ancestry with the African cluster. The nuclear C hybridization level calculated with IR formula was also low, with an average of 1.3%. These nuclear hybridization rates were much lower than the detected mitochondrial hybridization (average value 10.7%) ( Table 1).
To quantify the genetic differentiation, pairwise F ST values were obtained between La Palma population and A. mellifera subspecies. All the pairwise comparisons were statistically significant after Bonferroni correction. La Palma showed the lowest F ST values when compared to A. m. mellifera (F ST = .159 ± .043) and to A. m. iberiensis (F ST = .121 ± .025), although these two continental subspecies showed a closer genetic relationship between them (F ST = .028 ± .004). The comparisons of La Palma with African subspecies A. m. intermissa, A. m. major, and A. m. sahariensis showed high F ST values (F ST = .370 ± .102, F ST = .381 ± .118 and F ST = .401 ± .135, respectively), highlighting again the differentiation between Canarian and African bees, and the closer relationship between the Canarian honey bee and the M evolutionary branch. Finally, C and O lineage subspecies showed the highest F ST values when compared to La Palma (F ST = .572-.643).
To visualize the genetic relationship of La Palma honey bee with A. mellifera subspecies, a NJ tree was constructed based on D A distances for 6 microsatellites ( Figure 5(A)). South and Western Iberian populations were selected to build the tree. La Palma population grouped with Western European A. mellifera subspecies (M branch), clearly separated from African populations (A branch) and from C and O branches. To further explore the genetic relationship of Canarian populations with those of M branch, another NJ tree was built, based on data from 10 microsatellite loci from populations of La Palma, A. m. mellifera, and A. m. iberiensis ( Figure 5(B)). Canarian populations were grouped with Lisbon population with a high bootstrap value (99%). Aiming to evaluate the Iberian region as the origin of the La Palma honey bees, the A. m. iberiensis and the La Palma populations were also analyzed using the STRUCTURE software. All available populations from the Iberian Peninsula (Miguel et al., 2007) were included in the analysis. Results indicated the existence of two differentiated clusters, one of them with all the populations of the Iberian Peninsula and the other with those of La Palma (Figure 4(B) and Supplementary Figure 1  (B)). La Palma honey bees were located in a homogenous cluster, while the A. m. iberiensis populations from the southwestern region of the Iberian Peninsula (mainly, Lisbon) showed the lowest percentage of Iberian genetic component (Figure 4(B)). With regard to the classification of the individuals in the two clusters, 2.0% of La Palma honey bees were classified together with the Iberian ones, while 4.2% of the A. m. iberiensis (and specifically, 26.7% of those from Lisbon and 17.8% of those from Evora) were classified together with those from La Palma.

Geometric morphometrics
Since the inclusion of centroid size in the analyses did not modify the results and interpretations, all the following results were done using shape parameters alone, i.e., the tangent space projection matrix. Projections onto the first plane of the CVA performed at the population level (Figure 6) illustrated the large morphological gap between Canarian and African honey bees (CVA axis 1, 37.4%) and the close relationships between Canarian, Iberian and, more generally, M lineage honey bees, that were almost undistinguishable on the CVA axis 2 (23.0%). Both Canarian populations were largely mixed with the M populations, but occupied a somewhat marginal position on the second axis. Mahalanobis distances confirmed that the two most similar populations were those from Canary Islands (La Palma and Tenerife, D2 = 3.87). These populations showed greater distance values with M branch populations (with A. m. iberiensis, D2 = 6.24; with A. m. mellifera, D2 = 12.98) than those observed between A. m. mellifera and A. m. iberiensis (D2 = 4.70). Distances between African populations and both Western European and Canarian honey bees were higher, ranging from 22.02 to 27.56. Finally, looCV rates of colonies from La Palma showed 12.8% misclassified individuals, 11.4% attributed to A. m. iberiensis (one third of which corresponded to the population from Lisbon), 1.4% attributed to A. m. mellifera, and 0% to African subspecies. For the Tenerife population, the 12.0% misclassified individuals were all attributed to A. m. iberiensis. These high looCV correct classification percentages were only slightly lower than those observed in the A lineage populations (96%), while looCV percentages within the A. m. iberiensis subspecies ranged  from 0% to a maximum of 52%. All these results demonstrated that while Canarian honey bee populations undoubtedly belong to the A. m. iberiensis subspecies, they still exhibited a very substantial morphological differentiation from continental honey bees. Morphologically, the A. m. major and A. m. intermissa wings differed from the M and Canarian ones ( Figure 6, CVA axis 1), being broader with a large basal shift of the proximal wing part. The differences were somewhat less marked on the medial wing cells, but large shifts and displacements might be seen on the cubital and submarginal ones. The singularity of the Canarian honey bee was again confirmed by those marked wing pattern differences.

Discussion
The aim of this study was to inquire into the evolutionary history and genetic diversity of the Canarian honey bee in an attempt to provide useful information to assist in the management of the Canarian honey bee Conservation Project. To achieve this objective, a study of microsatellite markers, mitochondrial DNA, and geometric morphometric data was performed in a sample of 499 honey bee colonies of La Palma (the largest sampling to date).
Similarities on mtDNA haplotype frequencies have led De la Rúa et al. (1998Rúa et al. ( , 2001Rúa et al. ( , 2006 to suggest an African origin for the Canarian honey bee. However, the results obtained in the present study for mtDNA, microsatellite and geometric morphometrics were not in agreement with this conclusion, suggesting instead an Iberian origin for the Canarian honey bee. Geometric morphometrics reflected its close relationship with A. m. iberiensis, since no Canarian honey bees were classified in the A branch. With regard to microsatellites, F ST values, STRUCTURE analysis (Figure 4(A)) and NJ tree ( Figure 5(A)) also marked a clear distinction between La Palma honey bee and subspecies of the A, C, and O branches, given that all individuals from La Palma grouped with Western European A. mellifera subspecies (M branch). Mitochondrial data also pointed toward the same conclusion: COI-COII haplotype distribution of populations from La Palma appeared as a subset of Iberian haplotype distribution, since all haplotypes were also present in the Iberian Peninsula (Cánovas, De La Rúa, Serrano, & Galián, 2008;Garnery et al., 1998a;Miguel et al., 2007;Pinto et al., 2012Pinto et al., , 2013, while only three (A1, A8, A9) of the eight A haplotypes have been found in northwestern Africa (De La Rúa, Radloff, Hepburn, & Serrano, 2007;Franck et al., 2001). Moreover, the phylogenetic tree based on Cytb sequences ( Figure 2) showed a compact group formed by Iberian and Canarian colonies, and the median-joining network ( Figure 3) located all A1 haplotypes from African populations distant from Iberian and Canarian A1 haplotypes. These results indicated homoplasy for the A1 haplotype and suggested that the A I sublineage (Franck et al., 2001), which includes A1, A2, and A3 haplotypes, is not monophyletic. On the contrary, haplotypes included in A II and A III sublineages grouped closer than those of A I , and, although more research is necessary to verify the evolutionary pattern of COI-COII haplotypes, our results did not reject monophyly for these sublineages. In conclusion considering plausible homoplasy, the A I sublineage haplotype frequencies may not be adequate, per se, to infer the evolutionary history of the Canarian honey bee.
With respect to the exact location of the Canarian honey bee ancestors in the Iberian Peninsula, both nuclear and mitochondrial data suggested a Portuguese origin. Microsatellites showed a close phylogenetic relationship between Canarian honey bees and those from Portugal (Figures 4(B) and 5(B)). COI-COII haplotypes of mtDNA were also consistent with this relationship: in the Canary Islands the most common sublineage, is A III , which is predominant in the northern Portuguese populations . Moreover, Pinto et al. (2012) reported a high amount of novel A lineage haplotypes in Portugal, nearly all of them being from the A III sublineage, which would make Portugal the largest reservoir of A III haplotypes, and reinforce the idea of a Portuguese origin of the Canarian honey bee. It is worth pointing out that A III haplotypes are also frequent in the Portuguese archipelagos of Madeira and Azores (De La Rúa et al., 2006;, although they show different features regarding COI-COII haplotype distribution. The A14 haplotype (A III sublineage) is predominant in Azores, where no A1 (A I sublineage) is found, while A15 (A III sublineage) and in lesser degree A1 are the main haplotypes in Madeira (De La Rúa et al., 2006;Muñoz et al., 2013). The pattern found in La Palma is similar to the latter. With the caution that frequencies of only one marker are being considered, these results could indicate a common process of colonization of Madeira and Canary Islands. Indeed, Madeira is nearly midway between the Iberian Peninsula and the Canary Islands, while the Azores archipelago is at about 1000 km northwest of Madeira. Similar heterozygosity values in La Palma and Madeira (De La Rúa et al., 2006) do not clearly indicate which archipelago was first colonized, but the findings suggest a shared evolutionary history among these honey bees, similar to what was observed in the islands of the central Mediterranean region (Sheppard, Arias, Grech, & Meixner, 1997).
Historical and biological evidence shows that a major presence of African honey bees in the Canary Islands in the past was likely. In the first century AD, Pliny the Elder described a large supply of honey in the Canary Islands (Natural History of Pliny the Elder: Book 6, paragraph 37), which points toward the notion of bees on the islands prior to the Spanish conquest in the fifteenth century. Other insect evolution, such as Bombus terrestris canariensis (Rasmont, Coppee, Michez, & De Meulemeester, 2008;Widmer, Schmid-Hempel, Estoup, & Scholl, 1998) and Drosophila subobscura (Afonso et al., 1990) support the idea of an African origin for these "Plinyean" bees. However, in the sixteenth century, increasing apicultural activity occurred (Méndez, 2003). This fact could have led to a massive importation of Iberian honey bees, probably from Portugal, and might have derived in a complete replacement of the hypothetical ancient African honey bee gene pool. Indeed, microsatellites did not detect any trace of African honey bees in La Palma (Figure 4(A)). Nevertheless, it is worth pointing out that given the prevailing winds and sea currents (Juan, Emerson, Oromo, & Hewitt, 2000), a natural and ancient colonization from the Iberian Peninsula cannot be ruled out.
Despite the close phylogenetic relationship detected by the three kinds of markers between Iberian and Canarian populations, nuclear DNA F ST values and geometric morphometric results showed greater differences between Canarian honey bees and A. m. iberiensis than those observed between A. m. mellifera and A. m. iberiensis. In this sense, geographical isolation of the insular bees favors genetic drift and results in a loss of genetic diversity. In addition, micro-evolutionary processes in the Canary Islands, such as founder effects due to multiple introductions, bottleneck processes linked to volcanic activity on the islands and to the widespread infestation of Acarapis woodi in the 1950s (unpublished data), could have resulted in significant differentiation of the Canarian honey bees from A. m. iberiensis subspecies. Natural and artificial selection and also climatological and floral differences between the Canary Islands and Iberia probably produced a Canarian honey bee with differential adaptations. In fact, these honey bees show a different behavior with regard to Iberian ones, such as reduced defensiveness and a high hygienic aptitude (Gilles Fert, personal communication), and an absence of hibernation, due to the climatic conditions of the Canary Islands. For all these reasons, results obtained in the present study point toward the definition in the Canary Islands of an adapted ecotype within the M branch honey bee, i.e., a variant of A. m. iberiensis.
With regard to modern apiculture activity, massive importations of A. m. carnica and A. m. ligustica subspecies (C evolutionary branch) endanger the conservation of local populations, leading to the reduction and fragmentation of many populations (De La Rúa et al., 1998, 2001Franck et al., 1998, Garnery et al., 1998bMeixner et al., 2010;Soland-Reckeweg et al., 2009). Given this fact, local conservation initiatives have emerged in several European countries, such as the one on La Palma to protect the Canarian honey bee. The analysis of genetic diversity and screening of the hybridization with C lineage provide theoretical and practical genetic information for better animal management. The genetic singularity and local adaptation of the Canarian populations represent a significant component of the overall genetic diversity within the M branch. This branch presents lower diversity than A and C evolutionary branches, probably as a result of continuous regressions and post-glacial recolonization events (Miguel et al., 2007). So, from the perspective of conservation, the contribution of the Canarian honey bee to the M branch diversity highlights the need to preserve this genetic pool. In addition, the low heterozygosity of the Canarian populations, a negative for adaptation potential, strengthens the need of a special protection.
When screening for hybridization, nuclear data reflected very low values (1.3%), compared to the higher mitochondrial introgression value (10.7%). These differences in hybridization values might reflect a lower fitness of the foreign drones compared with the local drones and, as suggested by Muñoz, Pinto, and De La Rúa (2014), a partial reproductive isolation between native and introduced C subspecies. Detection of C haplotypes contrasts with the absence of mtDNA hybridization reported in La Palma a few years ago (De La Rúa et al., 1998, 2001, and could indicate a recent importation of queens to La Palma Muñoz et al., 2013). However, the present study provides valuable data for the Conservation Program, since hybridization appears to be restricted to a few populations in the middle and western part of the island. In addition, the hybridization level observed in La Palma remains low compared to other Canary Islands (De La Rúa et al., 1998, 2001Muñoz et al., 2013), and Muñoz et al. (2014) did not detect a temporal change at the nuclear level for . Soland-Reckeweg et al. (2009 detected hybrids in so-called "pure breeding areas" in Western European populations. They highlighted that conventional methods do not necessarily provide an effective barrier to gene flow between neighboring populations and suggested that management strategies should be adjusted for improved conservation management. In this sense, European legislation that is aimed at improving the conditions in which agricultural products are produced and marketed in Canary Islands (Council Regulation -EC, No 1454/2001), as well as the local law that forbids honey bee imports to La Palma and requires replacing foreign honey bees by local ones (Government of Canary Islands, 2001: BOC 49 of 20/4/2001, are necessary measures to prevent hybridization. Yet the choice of the North of the island as an Area of Special Protection seems appropriate to strengthen the barrier to hybridization. The orographic isolation of the populations of San Andrés and Barlovento turn it into an adequate refuge for the local genetic pool. Nevertheless, these two localities are not in Hardy-Weinberg equilibrium. Is it possible that the procedure of local drone saturation of this area, as a strategy to minimize the effect of imported queens, led to the detected lack of heterozygotes. To maximize genetic diversity in this area of Special Protection of the Local Honey Bee Conservation Program, the whole island seems to be an appropriate source of variability, albeit previously, hybridization with C lineage should be ruled out in candidate colonies.

Acknowledgments
We thank Carmela García Castaños, who put us in contact with Canarian Institutions. This manuscript was improved thanks to the comments of the editor and two anonymous reviewers.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This research was supported by the project "Análisis y estudio de DNA mitocondrial en abejas de la isla de La Palma," funded by the Gobierno de Canarias (Consejería de Agricultura, Ganadería, Pesca y Alimentación). Irati Miguel's work was supported by a PhD grant by the Basque Government (Departamento de Educació n, Universidades e Investigació n).

Supplemental data
Supplemental data for this article can be accessed here http:// dx.doi.org/10.1080/00218839.2016.1180017.