Evaluation of Morphometric Characters of Honeybee ( Apis mellifera L . ) Populations in the Lake Chad Basin in Central Africa

Samples of workers of honeybee were collected from 42 colonies in 13 localities in the Lake Chad Basin (parts of Cameroon, Chad, and Nigeria) and analysed using classical morphometry. Measurements of 35 morphological characters of body size, colour and pilosity were taken from 10 workers per colony and the data subjected to one-way analysis of variance (ANOVA), principal components analysis (PCA), hierarchical cluster analysis (HCA), stepwise discriminant analysis (DA) and Pearson’s product-moment correlation analysis. A one-way ANOVA revealed that means of 21 of the morphometric characters differed significantly (p < 0.05) among sampled localities while means of the remaining 14 characters did not (p > 0.05). The bees formed one cluster in a PCA. However, scatter plots of altitude against principal component 1 of PCA (loaded with characters of body size) revealed an increase of size of the bees along the gradient of the Lake Chad Basin. The coefficient of determination (R2) indicated that 88% and 77% of the variation in size might be explained by altitude in the southeastern and southwestern parts of the basin, respectively. Additionally, there was a very highly significant strong positive relationship between principal component 1 and altitude (r(30) = 0.618, p < 0.0005). Similarly, HCA and DA classified the colonies into three morphoclusters whose distribution closely followed the altitude of the area.


Introduction
Honeybees are of considerable economic importance, producing products of commercial value, such as honey and wax, and pollinating crops and wild plants.For example, Calderone [1] estimated the value of honeybees, through pollination of crops, to be US$11.68 billion in the United States of America, for the year 2009.Similarly, the contribution of honeybees to the economy of the United Kingdom, through pollination, was estimated at £191.80 million [2].The combined production of honey by the top 20 producer-countries for 2011 was estimated at 1.26 million metric tonnes valued at US$3.16 billion [3].In addition, products of honeybees, notably honey, royal jelly, propolis and bee venom, contribute to our well-being through their nutritional and therapeutic properties.
The natural range of the western honeybee, A. mellifera, is western Asia, Africa and Europe: From southern Scandinavia in the north to the Cape of Good Hope in the south, from Dakar in the west to the Urals, Mashhad and to the coast of Oman in the East.Geographical isolation and ecological adaptations resulted in the evolution of local populations showing considerable geographical variation, resulting in adaptation to local factors of climate, vegetation, pests and pathogens [4] [5].These adaptations may be lost due to human activities in beekeeping that affect wild honeybees in different ways: competition for floral sources, introduction of exotic genes, pests, parasites and diseases [6] [7].Therefore adequate knowledge of the natural diversity of local subspecies and ecotypes is essential for their management and conservation.To protect the biological diversity of local populations of honeybees in their natural habitats, these populations must, first of all, be characterised.One of the standard methods of characterising honeybees is classical morphometry [4].This method uses numeric data resulting from exact measurements of 36 morphological characters of body size, colour and pilosity from which means of colony characters are obtained for statistical analyses [8].Using this method, Ruttner [8] classified A. mellifera into 24 subspecies and four evolutionary lineages.However, whereas the European subspecies have been thoroughly studied, the study of their Asian and African counterparts is still in its infancy, in many places [4].With only 190 colonies, from 91 localities morphometrically analyzed [9], the western part of Africa (Countries in West and Central Africa, from Mauritania and Senegal in the west, to Chad in the east, then south to Namibia, through Zambia) is evidently under-studied.Although a few studies have been carried out, recently, in Nigeria [10]- [12], they have done little in improving the situation, due to their inadequacy in coverage and/or methodology.Thus, the present study attempted to improve our knowledge of the diversity of the honeybees of this region by analysing 42 colonies from 13 localities in the Lake Chad basin-through classical morphometry.
The main purpose of the study was to improve the availability of reference data in the region for future studies.

Description of the Area of Study
The area of study (Figure 1) lies approximately within 8˚ to 22˚ E and 8˚ to 14˚ N and covers parts of Chad, Cameroon and Nigeria.A summary of the important physical features of the area, based on Hepburn and Radloff [9] is given below.Most of the area lies in the Guinea and Sudan savanna zones though the lake itself lies in the Sahel, a transitional zone between savanna and desert.
The savannas consist of a mixture of grasslands and woody vegetation; moderate rainfall with a short rainy season (less than half of the year); and a high variation of mean annual temperature (10˚C to 15˚C).The density of vegetation, species richness, amount of rainfall and length of the rainy season decrease with the increase in latitude.Annuals flower at the end of the rainy season while trees flower during the dry season.
The Sahel is characterized by scanty rainfall; very short rainy season; very sparse vegetation of short trees and grass; and frequent droughts.
Altitude varies from 280 metres above sea level on the shores of the lake to about 1300 metres on the Jos Plateau.

Collection of Honeybees
Samples of workers of honeybee were collected from 42 colonies in 13 localities in Chad, Cameroon and Nigeria (Table 1 and Figure 1).Bees were collected, from wild nests or unmanaged traditional or top-bar hives populated by wild swarms, and preserved in 70% ethanol.The collection of samples and the morphometric data generated therefrom are stored at the Institut für Bienenkunde in Oberursel (Polytechnische Gesellschaft), University of Frankfurt, Germany.

Morphometric Measurements
Morphometric measurements of 35 characters were taken from 10 bees from each colony according to Ruttner [8] and Ruttner, Tassencourt and Louveaux [5].Measurements of hair and pigmentation were taken under a Leica dissecting microscope, fitted with an eyepiece graticle, at a magnification of 40×.Measurements of wings, legs and sternites were taken with the help of a Leica CCD camera connected to a desktop computer, using the measuring program Bee Morphometric, Version 1.02 [13].The details of the variables measured are shown in Table 2.

Statistical Analyses
First, the mean and standard deviation of each of the 35 morphometric characters were calculated for every colony.Then the means for colonies were used to calculate the means and standard deviations for localities.The data were subjected to a one-way analysis of variance (ANOVA) to compare different localities.Tukey HSD tests, on the means of the 35 characters, were used to detect significant differences between the localities.The level of statistical significance chosen was p = 0.05.Table 1.Localities from which samples of honeybee were collected for morphometric analysis (See map in Figure 1).A PCA, using colony means of 35 morphometric characters was run in order to detect any possible clusters.The suitability of PCA was assessed, prior to the analysis, using correlation coefficients of the variables, Kaiser-Meyer-Olkin (KMO) measure of sampling adequacy and Bartlett's test of sphericity [14] [15].Any variable that did not have a correlation with at least one other variable where r ≥ 0.3 should be removed from the analysis.The KMO measure is used as an index of whether there are linear relationships between the variables.Its value can range from 0 to 1, with values above 0.6 suggested as a minimum requirement for sampling adequacy.Bartlett's test of sphericity tests the null hypothesis that there are no correlations between any of the variables and, therefore, the variables cannot be reduced to a smaller number of principal components.For a PCA to be feasible the null hypothesis must be rejected: The result of Bartlett's test is used to take this decision.
Then a hierarchical cluster analysis was carried out with the assumption that samples from each locality formed a distinct population.Thus, 13 populations were assumed and the analysis was carried out using the means of the 35 morphometric characters for the 13 localities in the study.
A stepwise discriminant analysis (DA) was run in order to confirm the groups predicted by cluster analysis and to determine the discriminant characters.The suitability of DA was determined through log determinants and Box's M test.In DA the basic assumption is that the variance-co-variance matrices are equivalent.For this assumption to hold the log determinants should be equal and the Box's M test should not be significant.The latter tests the null hypothesis that the covariance matrices do not differ between groups formed by the dependent variable.Wilk's lambda was used to test the discriminatory power of the discriminant functions while the significance of the distance between group centroids was tested by F-statistic.
A Pearson's product-moment correlation analysis was run to assess the relationship between latitude, longitude and altitude and principal component 1 extracted by PCA (representing morphological characters loaded on it).Prior to analysis, the data were assessed for linearity, normal distribution and presence of outliers [14] [15].Interpretation of the magnitude of Pearson's correlation coefficient, r, was based on Cohen [16], that is: 0.1 < |r| < 0.3: small/weak correlation; 0.3 < |r| < 0.5: medium/moderate correlation; and |r| > 0.5: large/strong correlation.
All statistical analyses were carried out with IBM ® SPSS ® Statistics Version 20 with additional material from Burns and Burns [14] and Anonymous [15].
A one-way ANOVA revealed that means of 21 of the morphometric characters (length of cover hair on tergite 5, width of the tomentum band on the side of tergite 4, width of the dark stripe between the tomentum and the posterior rim of the tergite, length of femur, length of tibia, length of metatarsus, width of metatarsus, longitudinal diameter of tergites 3 and 4, longitudinal diameter of sternite 3, transversal diameter of wax plate of sternite 3, transversal diameter of sternite 6, length of fore-wing, five angles of wing venation (D7, E19, J16, K19 and N23), pigmentation of tergite 3, pigmentation of scutellum and pigmentation of plates of scutellum) differed significantly (p < 0.05) between sampled localities while means of the remaining 14 characters (pigmentation of the second and fourth tergites, longitudinal diameter of wax plate (mirror) of sternite 3, distance between wax mirrors of sternite 3, longitudinal diameter of sternite 6, width of fore-wing, the two cubital distances, and angles A4, B4, G18, J10, L13, and 026) did not (p > 0.05).
In order to investigate the similarity of the honeybee colonies under study, a PCA, using colony means of six morphometric characters (29 characters were excluded from the analysis during a preliminary PCA due to their   failure to meet some conditions) of 10 worker honeybees from each of 42 colonies at 13 localities, was run to detect possible clusters.The suitability of PCA was assessed prior to analysis.Inspection of the correlation matrix showed that all variables had the minimum requirement of at least one correlation coefficient greater than 0.3.The overall Kaiser-Meyer-Olkin (KMO) measure was 0.68 with individual KMO measures from 0.53 to 0.82, thus meeting the minimum requirement for sampling adequacy.Bartlett's test of sphericity was statistically significant (p < 0.0005), suggesting that the data could be appropriately analysed using PCA [14] [15].Three principal components, with eigenvalues 2.97, 1.23 and 0.72 each and accounting for 81.57% of the total variance were extracted.A Varimax orthogonal rotation was employed to aid interpretability.There were strong loadings of characters of size on all three components and pigmentation of plates of scutellum on component 2 (Table 7).
As may be seen in Figure 2, scatter plots of various combinations of the three principal components did not produce distinct clusters.However, plots of altitude against principal component 1 revealed an increase of size of the bees along the gradient of the Lake Chad Basin.The coefficient of determination (R²) indicated that 88% and 77% of the variation in size might be explained by altitude, respectively, in the southeastern (Figure 3(a)) and southwestern (Figure 3(b)) parts of the basin.
Since PCA did not produce a clear clustering of the colonies, an alternative method, hierarchical cluster analysis, using means of 35 morphometric characters for the 13 localities under investigation, was used.Based on the output of this analysis (Figure 4(a) and Table 8), three tentative clusters were assumed: Cluster 1 made up of colonies from Bauchi, Jos, Maiadua, Mubi, Touboro and Koumra; cluster 2 colonies from Afunori, Maga and N'Djamena; and cluster 3 colonies from Am Timan, Bousso, Sarh, Moundou.
A stepwise discriminant analysis (DA) was run, in order to predict the membership of the 42 colonies of A. mellifera, among the three a priori groups defined by the cluster analysis, as a means of confirming these groups.35 morphometric characters of workers were used as predictor variables.Significant mean differences (ANOVA; p < 0.05) were observed for all, but 12, of the variables.While the log determinants were similar (−9, −10 and −11, respectively, for the three groups), Box's M indicated that the assumption of equality of covariance matrices was not violated (p > 0.05).Thus, DA was deemed appropriate [14] [15].
Two canonical discriminant functions (with eigenvalues 4.6 and 3.4, respectively) were used in the analysis.They explained 57.8% and 42.2% of the total variance, respectively, and their Wilk's Lambda values, assessed by      A. mellifera from the Lake Chad basin by hierarchical cluster analysis (using average linkage (between groups)) using means of 35 morphometric characters for the 13 localities.Based on the dendrogram, three tentative groups were identified for entry into stepwise discriminant analysis; (b) Confirmation of the three groups by discriminant analysis.95% of cross-validated grouped cases were correctly classified into their original groups and the distances between group centroids were highly significant (p < 0.0005) according to F-statistic.75% confidence ellipses are shown.chi-square, were highly significant (p < 0.0005), suggesting they had sufficient discriminatory power to group the cases.The discriminant variables used in the analysis and their correlations with the discriminant functions are shown in Table 9.
Overall, 95% of cross-validated grouped cases were correctly classified.13 of the 15 cases in cluster 1 were classified with a posterior probability of at least 95% while 12 of the 13 cases in cluster 2 were also classified  with the same posterior probability.For cluster 3 all, but one of the 13 cases, were as well classified with a posterior probability of 95% and above.
A pairwise comparison of the groups, using F-statistic, revealed a very highly significant difference between the groups' centroids (p < 0.0005).A scatter plot of the two discriminant functions is shown in Figure 4(b) and the distribution of the members of the three groups (morphoclusters), in the area of study, in Figure 5.
A Pearson's product-moment correlation analysis was run to assess the relationship between latitude, longitude and altitude and principal component 1 extracted by PCA and loaded with characters of body size.Preliminary analyses showed the relationship to be linear with all variables normally distributed (as assessed by visual inspection of normal Q-Q plots) and all outliers were removed [15].There was a very highly significant strong positive relationship between principal component 1 and altitude (r (30) = 0.618, p < 0.0005) and a significant moderate negative relationship between this component and longitude (r (30) = −0.391,p < 0.01).However, this component had a non-significant weak negative relationship with latitude (r (30) = −0.281,p > 0.05).
In order to isolate the effect of longitude on the relationship between altitude and principal component 1, and vice vasa, the relationship of either variable with principal component 1 was subjected to a first-order partial correlation, controlling for the effects of the other.The correlation between altitude and principal component 1, controlling for longitude, was statistically significant (r (29) = 0.371, p < 0.05) indicating that a relationship between altitude and principal component 1 exists above and beyond the effects of longitude.On the other hand, the correlation between longitude and this component, controlling for altitude, was not statistically significant (r (29) = 0.255, p > 0.05), indicating that a relationship between longitude and principal component 1 does not exist above and beyond the effects of altitude.This underscores the importance of altitude in the variation of size of bees in the Lake Chadic basin (Figure 3).

Discussion
As may be seen in Table 10, the mean values of a set of the morphometric characters measured in this study are in general agreement with those reported for the subspecies of A. mellifera in sub-Saharan Africa [8] [12], except some suspicious values reported by Ajao, Oladimeji [10] and Oyerinde, Dike [11].This agreement is a confirmation of the correctness of the measurements taken in this study.
The trend observed in this study, in which the size of the bees increases with altitude, was similarly reported from East Africa [17] [18] and it may be due to the effect of an ecocline -a graduation in measurable characters observed in zones of different altitudes [19] [20]-according to Bergmann's rule: larger animals of the same species are found at higher altitudes or latitudes.This observation is supported by the strong positive correlation between size and altitude reported in this study.

Conclusion
Based on the agreement between the results of this study and those of previous studies, it is concluded that the honeybees of this area are morphometrically pure populations of sub-Saharan A. mellifera.

Figure 1 .
Figure 1.The area of study (enclosed by the black thick line) showing sampled localities.It comprises southern part of Chad, northern part of Cameroon and northeastern part of Nigeria.Inset: A map of Africa showing the location of the area of study.

Figure 2 .
Figure 2. PCA plots using the colony means of six morphological characters of workers of A. mellifera from 13 localities in the Lake Chad basin.All three components were loaded with characters of body size.In addition component 2 was loaded with pigmentation of plates of scutellum.The colonies are coded according to their localities.

Figure 3 .
Figure 3. Variation of size of A. mellifera in the Lake Chad basin, illustrated by scores of the first principal component of PCA (loaded with characters of body size).Size increases with altitude.The coefficient of determination (R²) indicates that 88% and 77% of the variation in size may be explained by altitude, respectively, in the southeastern (A) and southwestern (B) parts of the basin.

Figure 4 .
Figure 4. (a) Classification of 13 populations of A. mellifera from the Lake Chad basin by hierarchical cluster analysis (using average linkage (between groups)) using means of 35 morphometric characters for the 13 localities.Based on the dendrogram, three tentative groups were identified for entry into stepwise discriminant analysis; (b) Confirmation of the three groups by discriminant analysis.95% of cross-validated grouped cases were correctly classified into their original groups and the distances between group centroids were highly significant (p < 0.0005) according to F-statistic.75% confidence ellipses are shown.

Figure 5 .Table 10 .
Figure 5. Distribution of members of three morphoclusters of A. mellifera from the Lake Chad basin defined by hierarchical cluster analysis and confirmed by stepwise discriminant analysis.Only colonies assigned with a posterior probability of at least 95% (90% of total number of colonies in the study) are shown.Each circle represents one colony.The black thick line demarcates the area of study.Inset: A map of Africa showing the location of the area of study.Table10.Comparison of values (mean ± s.d.) of some morphometric characters of subspecies of A. mellifera from sub-Saharan Africa from various sources.

Table 2 .
List of characters measured for morphometry1.

Table 3 .
Means and standard deviations of morphological characters of body hair (mm) and pigmentation of 10 worker bees each from (N) colonies from 12 localities in the Lake Chad basin.

Table 4 .
Means and standard deviations (mm) of characters of the hind leg of 10 worker bees each from (N) colonies from 12 localities in the Lake Chad basin.
*** Very highly significant (p < 0.001), according to one-way ANOVA.Means, within a column, followed by the same letter are not significantly different at 0.05 according to Tukey HSD test.Numbers in brackets are Ruttner numbers.

Table 5 .
Means and standard deviations (mm) of characters of the abdomen of 10 worker bees each from (N) colonies from 12 localities in the Lake Chad basin.

Table 6 .
Means and standard deviations of characters of the fore-wing of 10 worker bees each from (N) colonies from 12 localities in the Lake Chad basin.
* Significant (p < 0.05); ** highly significant (p < 0.01); *** very highly Significant (p < 0.001) and † not significant (p > 0.05), according to one-way ANOVA.Means, within a column, followed by the same letter are not significantly different at p = 0.05 according to Tukey HSD test.Measurements of distance are in mm and of angles in degrees.Fwl: Fore-wing, length; Fww: Fore-wing, width; Cub-a, b: Cubital vein, distance a, b; A4-O26: 11 angles of wing veins.Numbers in brackets are Ruttner numbers.

Table 7 .
Rotated component matrix for PCA, with Varimax rotation, of morphometric characters of A. mellifera colonies from the Lake Chad basin.Major loadings for each character are shown in bold.

Table 8 .
Proximity (dissimilarity) matrix of 13 populations of A. mellifera from the Lake Chad basin.

Table 9 .
Pooled within-groups correlations between discriminating variables and standardized canonical discriminant functions.The largest absolute correlation between each variable and any discriminant function is shown in bold.