Geometric morphometric discrimination of the three African honeybee subspecies Apis mellifera intermissa, A. m. sahariensis and A. m. capensis (Hymenoptera, Apidae): Fore wing and hind wing landmark configurations

Choukri Barour; Michel Baylac

doi:10.3897/jhr.52.8787

Abstract

In this study, a landmark-based geometric morphometric analysis was carried out on three honeybee subspecies: Apis m. intermissa and A. m. sahariensis collected from Algeria, and, as a reference, A. m. capensis collected from South Africa. The aim of this study was to discriminate honeybee subspecies by patterns of shape variation of fore and hind wings. A total of 540 wings from 270 honeybee workers were analyzed. Our results revealed very high cross-validation classification rates (96.7% based on fore wing shape and 99.6% based on the combination of fore and hind wing forms respectively). Discrimination was better using shape and form (shape + centroid size) of the fore wings than of the hind wings. The wing form parameters were found to differ significantly in shape and centroid size among the three analyzed subspecies. Finally, it may be concluded that landmark-based geometric morphometrics could be a powerful tool to characterize the Algerian honey bees.

Keywords

Bees, Africa, biodiversity, wing venation, morphometry, Procrustes superimposition, canonical variate analysis (CVA)

Introduction

For the African honeybees only limited studies are available. For example, in Algeria only a few studies were published on the local subspecies. Barour et al. (2005) by using distance measurements revealed a south-to-north clinal increase in overall size with significant population differences. This study filled a former gap in the knowledge of Apis mellifera intermissa population variability between Morocco and Tunisia. Recently, Barour (2012), Barour et al. (2011), Chahbar et al. (2013) and Loucif-Ayad et al. (2015) have studied the diversity of the local Algerian honeybees by using morphometric and molecular data. On the one hand, Barour et al. (2011) have demonstrated that populations can be distinguished on the basis of the shape of the fore wing venation of A. m. intermissa workers. On the other hand, Chahbar et al. (2013) have shown that naturally distributed honeybees in Algeria correspond to the African evolutionary lineage by determining their mitochondrial haplotype and the variation of ten microsatellite loci. Loucif-Ayad et al. (2015) have assessed the phylogenetic and population structure of many populations and have confirmed the African origin of these A. m. intermissa and A. m. sahariensis Algerian populations.

In this study, we wish to take part in an ongoing inventory program for setting up a morphological database of the Algerian honey bees. The purpose was to distinguish between African honeybee specimens: two local subspecies from Algeria, A. m. intermissa (Buttel-Reepen, 1906), A. m. sahariensis (Baldensperger, 1924), and one for comparison, A. m. capensis (Eschscholtz, 1821) from South Africa. Landmarks were recorded on fore and hind wing and analyzed using generalized least squares Procrustes analysis (GPA) to examine the discrimination between subspecies. Unfortunately, our A. m. sahariensis and A. m. capensis sample sizes are relatively small because in Algeria the populations of A. m. sahariensis are declining in abundance. They are under conservation plans and their geographic distribution is very restricted. A sample of A. m. capensis was provided by our South African colleagues and comprised only 90 bees. For these reasons, we have adapted the A. m. intermissa sample size to those of A. m. sahariensis and A. m. capensis, and we have not taken into account the biogeographic information.

Material and methods

Sample collection

Samples of workers of A. m. intermissa and A. m. sahariensis from Algeria and A. m. capensis from South Africa were collected from fixed-site beehives (no transhumance activities) and preserved in alcohol. The honeybees of A. m. sahariensis were provided by professional beekeepers and collected in the same eco-climatic region (almost a Saharan climate). The honeybees of A. m. intermissa were collected in the same eco-climatic region (northeastern) from two adjacent apiaries. Our A. m. intermissa samples were identified by a semi-automatic expert identification system of bee workers (Baylac et al. 2008) as A. m. intermissa, which is the only abundant honeybee subspecies in the northern regions of Algeria. By contrast, A. mellifera sahariensis ranges from Aïn Sefra and Béchar in Algeria through the oases of the Sahara south of the Atlas Mountains, its natural distribution is thus restricted to few Algerian Saharan regions. Incidentally, there are no breeding activities of A. m. sahariensis in the northern regions of Algeria. As a result, the collected A. m. intermissa samples (northeastern in our case) are very distant (approximately 850 km) from those of A. m. sahariensis (closest to the Sahara). Random samples of 15 workers per colony were used to build the final datasets that comprise 270 worker honeybees, i.e., a total of 540 fore and hind wings from the three subspecies. Geographic data of our samples are mentioned in Table 1.

Table 1.

Download as

CSV

XLSX

Origin and sample size for Apis mellifera subspecies.

Subspecies	Colony no.	Bee no.	Locality and province	Geographic coordinates
A. m. intermissa	03	45	Guelma-Aïn Makhlouf (Algeria)	36°14'28"N 7°14'39"E
A. m. intermissa	03	45	Guelma-Aïn Arbi (Algeria)	36°15'53"N 7°23'59"E
A. m. sahariensis	03	45	Nâama-Aïn Sefra (Algeria)	32°44'33"N 0°34'50"W
A. m. sahariensis	03	45	Béchar (Algeria)	31°36'58"N 2°13'56"W
A. m. capensis	06	90	Grahamstown (South Africa)	33°18'57"S 26°31'10"E
Total	18	270 (540 wings)

Data acquisition

For data acquisition, the right fore wing (FW) and hind wing (HW) of workers were cut at their base. Then the wings were temporarily slide mounted in distilled water and photographed using a digital camera attached to a Z6 APO Leica macroscope. For the FW, the coordinates of nineteen homologous landmarks (LM) were recorded as defined by Smith et al. (1997) (Fig. 1). Seven LM of the HW were then defined; the 1^st and the 7^thLM were defined as the result of the intersection of the first and the last hamuli with the radial vein (Fig. 1). All LM coordinates were digitized twice using the TPSDig2 software version 2.10 (Rohlf 2008) by the first author. The average of both repetitions was used in the final analyses in order to minimize the measurement error. All the LM correspond to type I LMsensuBookstein (1991) except for LM15 on FW (maximum curvature of veins, a type III LM). The same set or a slightly reduced set of LMs has been employed in most GPA studies of honeybees (Baylac et al. 2008, Miguel et al. 2011, Barour et al. 2011; Kandemir et al. 2011, Oleksa and Tofilski 2015), but other authors have also used HW shape data from the A. m. intermissa and A. m. sahariensis populations (Barour 2012) and from A. mellifera subspecies (Dolati et al. 2013).

Figure 1.

Location of the landmarks digitized on a right fore and hind wing of Apis mellifera workers (drawn to the same scale). MR: marginal cell, CC: cubital cell, MC: median cell, SMC: sub-median cell, and RC: radial cell.

Data analysis

Differences in FW and HW shapes between the three studied subspecies were investigated by the `Rmorph` library (Baylac 2007). The LM coordinates were then aligned using the Generalized least squares Procrustes Analysis (GPA) algorithm (Goodall 1995). The FW and HW sizes were estimated by the log-transformed centroid size (CS) (Bookstein 1991). Differences between the three honeybee subspecies were analyzed by canonical variate analyses (CVA). Misclassification rates for both wing shapes and forms (shape + CS) were estimated using a leave-one-out cross validation procedure (6 contrasts were used from shape and form data and their mixing between FW and HW). Wing size differences among honeybee subspecies and colonies were analyzed through box plots and one- and two-way ANOVAs of CS. Pairwise comparisons used t-tests with Holm’s p-value adjustment (Rice 1989). Statistical analyses were performed using R, version 3.0.3 (R Core Team 2014; Ihaka and Gentleman 1996) for Windows (http://cran.r-project.org/).

Results

Centroid size differences between the three A. mellifera subspecies

The analysis of size is presented for the logarithm of centroid size (CS). For the FW, a two-way ANOVA revealed highly significant differences between the honeybees for both the subspecies (p < 0.001) and the colony (p < 0.001) levels. In the case of the HW, the differences between the same honeybee samples were also highly significant (p < 0.001) for subspecies and colony levels. Furthermore, the subspecies × colony interaction term of the two-way ANOVA was highly significant (p < 0.001) for both the FW and HW. Pairwise comparisons using t-tests with Holm’s p-value adjustment showed that CS_FW differed significantly between intermissa-capensis (p < 0.001), sahariensis-capensis (p < 0.001) and intermissa-sahariensis (p < 0.001). Pairwise comparisons were also conducted for CS_HW and the subspecies were found to be statistically different only between intermissa-capensis (p < 0.001) and sahariensis-capensis (p < 0.001), but no significant difference was established between the Algerian subspecies intermissa-sahariensis (p > 0.05) (Fig. 2).

Figure 2.

Box plots of the logarithm of wing centroid size for each subspecies of Apis mellifera.

Wing shape and form differences

For both FW and HW, the MANOVA results of the shape data were highly significant (p < 0.001) between subspecies and also between colonies. The interaction term between the two factors (levels) subspecies × colony was also highly significant (p < 0.001) for both Shape_FW and Shape_HW, which revealed that more complex patterns existed within each honeybee subspecies.

Concerning Shape_FW variability, a scatterplot of the first two canonical variates (CV) (Fig. 3A) indicated a good discrimination of subspecies. Apis m. intermissa, and A. m. capensis were entirely separated by the first CV (explaining 57.95% of the variance), while A. m. sahariensis was largely separated along the second CV (explaining 42.05%). In contrast, for Shape_HW the three subspecies were largely overlapping along the first and the second CV (explaining 61.56% and 38.44%, respectively; Fig. 3B).

Figure 3.

Shape variability among Apis mellifera intermissa, A. m. sahariensis and A. m. capensis: first two canonical variates. A fore wing shape B hind wing shape.

FW and HW shape deformations along the first two CVs are visualized in Figures 4A, B, C and D. The first CV (Fig. 4A) mainly described the extreme shape differences observed between A. m. intermissa and the two others subspecies, notably A. m. capensis. Quite noticeably, the FW shape deformations along this CV involved almost the whole set of wing cells. For example, in A. m. intermissa (Black line in Fig. 4A) the marginal cell appeared longer and narrower than in A. m. capensis (Grey line in Fig. 4A). The first cubital cell also appeared narrower in A. m. intermissa than in A. m. capensis. By comparison, along the first CV the A. m. sahariensisFW shape patterns are more similar to those recorded in A. m. capensis. Shape changes associated with the second CV (Fig. 4B) correspond to differences between A. m. sahariensis (Black line in Fig. 4B) and the two others honeybee subspecies (Grey line in Fig. 4B). In addition, also on the second CV, the FW width of A. m. intermissa and A. m. capensis appear very similar and less wide than in A. m. sahariensis. Figure (4C) demonstrates the HW shape differences and the associated deformations of slightly overlapping clusters. On the whole, HW shape deformations along the first CV (Fig. 4C) was affected by the shifting of all seven LMs with contraction of both the radial cell and the radial vein (Distance LM1-LM7) in A. m. intermissa (Grey line in Fig. 4C). In comparison with A. m. intermissa, A. m. sahariensis is characterized by a shorter median vein (Black line in Fig. 4C) and a longer radial vein. Moreover, the HW shape of A. m. capensis has an intermediate LM configuration. In general, the contraction of both distances LM1-LM2 and LM2-LM7 could be interpreted as a variation in the number of hamuli among the three honeybee subspecies. Finally, the deformations along the second CV (Fig. 4D) highlight the major HW shape differences and LM shifts between the two Algerian subspecies (Grey line in Fig. 4D) and A. m. capensis (Black line in Fig. 4D).

Figure 4.

Extreme shape differences between Apis mellifera intermissa, A. m. sahariensis and A. m. capensis along the first two canonical variates (Fig. 3A, B). A and B fore wing shape differences along the first and second canonical variate, respectively. C and D hind wing shape differences along the first and second canonical variate, respectively (scale factor ×3 and ×2 respectively). Grey lines depict the shape associated with the negative values and black lines the shape associated with the positive values of the respective canonical variate.

Moreover, Table 2 shows the classification rates in the cross-validation tests on wing shapes and forms (mixing the logarithm of CS with shape parameters) among the A. m. intermissa, A. m. sahariensis and A. m. capensis. The cross-validated classifications calculated on Shape_FW and Shape_HW parameters assigned correctly 96.66% and 77.40% of the bees respectively; while cross-validated classifications on form (Shape + CS) reached 99.63% for FW and 90.00% for HW. Two other combination types between the FW and HW were used to verify that these combinations could enable even clearer discrimination. First, a high rate of correct classification 99.63% was obtained by using the FW and HW form combination. Second, we observed another high rate of discrimination of 97.41% by using a combination of FW and HW shape.

Table 2.

Download as

CSV

XLSX

Reclassifications results at the subspecies level: A. m. intermissa, A. m. sahariensis, and A. m. capensis (always n = 90).

Type of wing	Variables	Reclassification rate
FW	Shape_FW	96.66%
HW	Shape_HW	77.40%
FW	Shape_FW + CS_FW	99.63%
HW	Shape_HW + CS_HW	90.00%
FW + HW	Shape_FW + Shape_HW	97.41%
FW + HW	Shape_FW + CS_FW + Shape_HW + CS_HW	99.63%

Discussion

Morphometric methods have been advantageously applied to evaluate biodiversity and for taxonomic purposes (Francoy et al. 2011). In the present study, the exploratory analyses on FW and HW shape of African honeybees yielded three well-defined clusters with high reclassification percentages (96.66% for FW shape and 99.63% for a combination of FW and HW forms). The results indicated that geometric morphometrics using landmarks efficiently distinguished A. m. intermissa, A. m. sahariensis and A. m. capensis. This result is supported by several other publications on A. mellifera (Baylac et al. 2008, Tofilski 2008, Barour et al. 2011, Francoy et al. 2011, Miguel et al. 2011). Kandemir et al. (2011) also reported that a landmark analysis of wing shape could be used as a reliable tool to discriminate among honeybee subspecies. Furthermore, Oleksa and Tofilski (2015) reported that in some studies, morphometrics proved to be even more effective in the identification of subspecies than molecular markers, and that the morphological characters were also more suitable for distinguishing ecotypes within A. mellifera subspecies. In contrast, Dolati et al. (2013) reported that for A. m. meda colonies of many Iranian populations only 68.2% were correctly classified by using the FW shape and 43% by using the HW shape.

Our analysis of the HW venation also showed an important shift in the position of the 1^st, 2^nd and 7^th landmark. This variability could be explained by the distribution of the number of hamuli within A. mellifera. Indeed, the number of hamuli can be promising for distinguishing subspecies of A. mellifera, as the number of hamuli and their extent on the edge of the HW of honeybees have high heritability values and are readily modified by genetic selection (Hepburn and Radloff 2004).

In conclusion, the results presented here showed that (i) A. m. intermissa, A. m. sahariensis and A. m. capensis populations could be distinguished on the basis of the shape of worker wing venation. There was only a slight increase in the classification rate with form parameters, indicating that the differences involved mostly shape, whereas size was rather less important; (ii) FW venation was more powerful for discrimination than HW venation.

Acknowledgements

The authors thank L. Deharveng (Directeur UMR CNRS 7205), R. Cornette and all the people that work in the Morphometrics platform, MNHN, Paris, France. They thank the editors of JHR and the three reviewers for their comments that improved the final version of this paper. They are also grateful to M. Hamzaoui and many beekeepers for having facilitated the collection of the Algerian honeybee colonies. In addition, they profoundly thank Pr. R. Hepburn, Pr. S. Radloff and Dr. V. Dietemann for the Cape honeybee samples. They are especially grateful to Pr. S. Radloff (Rhodes University) for the linguistic corrections to the manuscript. This research work was granted by an Exceptional National Program fellowship from the Algerian Ministry of Higher Education and Scientific Research.

References

Barour C (2012) Analyse de la Biodiversité des Populations d’Abeilles Mellifères Apis mellifera intermissa (Buttel-Reepen, 1906) (Hymenoptera: Apidae) dans le Nord Algérien: Morphométrie Moderne Basée sur la Configuration des Points-Repères (Landmarks). PhD Thesis, University of Badji Mokhtar, Annaba.

Barour C, Tahar A, Radloff SE, Hepburn HR (2005) Multivariate analysis of honeybees, Apis mellifera Linnaeus (Hymenoptera: Apidae) of the northeastern and southern regions of Algeria. African Entomology 13(1): 17–23.

Barour C, Tahar A, Baylac M (2011) Forewing shape variation in Algerian honeybee populations of Apis mellifera intermissa (Buttel-Reepen, 1906) (Hymenoptera: Apidae): A landmark based geometric morphometrics analysis. African Entomology 19(1): 11–22. doi: 10.4001/003.019.0101

Baylac M (2007) Rmorph: a R geometric and multivariate morphometrics library. [Available from the author at baylac@mnhn.fr]

Baylac M, Garnery L, Tharavy D, Pedraza-Acosta J, Rortais A, Arnold G (2008) ApiClass, an automatic online wing morphometric expert system for honeybee worker identification. http://apiclass.mnhn.fr [accessed on 22 October 15]

Bookstein FL (1991) Morphometric Tools for Landmark Data: Geometry and Biology. Cambridge University Press, Cambridge, 456 pp.

Chahbar N, Muñoz I, Dall’Olio R, De la Rúa P, Serrano J, Sallaheddine D (2013) Population structure of North African honey bees is influenced by both biological and anthropogenic factors. Journal of Insect Conservation 17: 385–392. doi: 10.1007/s10841-012-9520-1

Dolati L, Rafi JN, Khalesro H (2013) Landmark-based morphometric study in the fore and hind wings of an Iranian race of European honeybee (Apis mellifera meda). Journal of Apicultural Science 57(2): 187–196. doi: 10.2478/jas-2013-0028

Francoy TM, Grassi ML, Imperatriz-Fonseca VL, May-Itzá W, Quezada-Euán JJG (2011) Geometric morphometrics of the wing as a tool for assigning genetic lineages and geographic origin to Melipona beecheii (Hymenoptera: Meliponini). Apidologie 42(4): 499–507. doi: 10.1007/s13592-011-0013-0

Goodall CR (1995) Procrustes methods in the statistical analysis of shape revisited. In: Mardia KV, Gill CA (Eds) Current issues in statistical shape analysis. University of Leeds Press, Leeds, 18–33.

Hepburn HR, Radloff SE (2004) The wing coupling apparatus and the morphometric analysis of honeybee populations. South African Journal of Science 100: 565–570.

Ihaka R, Gentleman R (1996) R: a language for data analysis and graphics. Journal of Computational and Graphical Statistics 5: 299–314.

Kandemir I, Özkan A, Fuchs S (2011) Reevaluation of honeybee (Apis mellifera) microtaxonomy: a geometric morphometric approach. Apidologie 42: 618–627. doi: 10.1007/s13592-011-0063-3

Loucif-Ayad W, Achou M, Legout H, Alburaki M, Garnery L (2015) Genetic assessment of Algerian honeybee populations by microsatellite markers. Apidologie 46: 392–402. doi: 10.1007/s13592-014-0331-0

Miguel I, Baylac M, Iriondo M, Manzano C, Garney L, Estonba A (2011) Both geometric morphometrics and microsatellite data support the differenciation of the Apis mellifera M evolutionary branch. Apidologie 42(2): 150–161. doi: 10.1051/apido/2010048

Oleksa A, Tofilski A (2015) Wing geometric morphometrics and microsatellite analysis provide similar discrimination of honey bee subspecies. Apidologie 46: 49–60. doi: 10.1007/s13592-014-0300-7

R Core Team (2014) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.r-project.org/

Rice WR (1989) Analyzing tables of statistical tests. Evolution 43: 223–225. doi: 10.2307/2409177

Rohlf F (2008) tpsDig2, version 2.10. Department of Ecology and Evolution, State University of New York, Stony Brook. http://life.bio.sunysb.edu/morph [accessed on 22 October 15]

Smith DR, Crespi BJ, Bookstein FL (1997) Fluctuating asymmetry in the honey bee, Apis mellifera: effects of ploidy and hybridization. Journal of Evolutionary Biology 10: 551–574. doi: 10.1007/s000360050041

Tofilski A (2008) Using geometric morphometrics and standard morphometry to discriminate three honeybee subspecies. Apidologie 39: 558–563. doi: 10.1051/apido:2008037