Geometric morphometric discrimination of the three African honeybee subspecies Apis mellifera intermissa

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.


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.

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.

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 th LM 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 LM sensu Bookstein (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 andTofilski 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).

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).

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).

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).
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. sahariensis FW 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).
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.

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.ers 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.
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/).

Figure 1 .
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.

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

Figure 3 .
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.

Figure 4 .
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.

Table 1 .
Origin and sample size for Apis mellifera subspecies.