Computer software for identification of honey bee subspecies and evolutionary lineages

Within the western honey bee (Apis mellifera), there are more than 20 recognised subspecies. It is well known that these subspecies differ in their wing venation patterns. However, there is a demand for efficient tools to identify honey bee subspecies, ecotypes, populations or hybrids. The aim of this study was to develop a fast and easy identification method based on analysing forewing vein patterns of honey bees by geometric morphometrics. Reference samples for the subspecies were obtained from the Morphometric Bee Data Bank in Oberursel, Germany. These contained 187 honey bee colonies allocated into 25 subspecies from four evolutionary lineages. The identification of evolutionary lineages of honey bees based on forewing venations proved to be highly reliable, which confirms earlier studies. The accuracy of honey bee subspecies identification was less consistent and ranged from 100 to 50% and was particularly low in African honey bees. The obtained identification data were exported to the IdentiFly computer software, which is freely available.


INTRODUCTION
The western honey bee (Apis mellifera ) has a wide geographical range, which covers almost the entire area of Europe, Africa, the Near East and central Asia (Ruttner 1988;Sheppard and Meixner 2003;Chen et al. 2016). Within this range, there are markedly different climates and environments. Such diverse conditions and the history of spread and isolation of subpopulations resulted in a notable variation of morphological traits. These traits, however, are rarely distinct in single subspecies, and multivariate statistical approaches were indispensable for reliable identification. By incorporating behavioural traits, numerous subspecies were defined (last reviewed by Engel 1999). These subspecies are grouped into four evolutionary lineages based on analysis of morphological characteristics (Ruttner 1988). Lineage A occurs in central and southern Africa, lineage C in southwest Europe, lineage M in Northern and Western Europe and lineage O in the Middle East (Kandemir et al. 2011;Meixner et al. 2013). Later molecular analyses mostly confirmed the morphological groupings, with few minor discrepancies (Franck et al. 1998;Alburaki et al. 2011).
Discrimination of the subspecies is important for the conservation of the honey bee biodiversity. Because of honey bee queen trade and migratory beekeeping, the natural ranges of subspecies are increasingly disturbed (De la Rúa et al. 2009). Furthermore, feral honey bee populations shrank Electronic supplementary material The online version of this article (doi:10.1007/s13592-017-0538-y) contains supplementary material, which is available to authorized users. or vanished in certain parts of Europe (Moritz et al. 2007). Monitoring local subspecies and limiting the introduction of non-native subspecies would help to conserve honey bee genetic biodiversity (De la Rúa et al. 2009). Moreover, discrimination of honey bee subspecies can be important in zoogeographic studies. In some parts of Asia and Africa, the distribution ranges of honey bee subspecies have not been fully explored (Meixner et al. 2013). Precise zoogeographical maps based on standardised identification methods would make an important contribution to future conservation programs.
Proper identification of the subspecies and evolutionary lineages is also important in honey bee breeding. Breeders often declare that their breeding lines belong to a particular subspecies. In order to keep those breeding lines as pure stocks, breeders need to eliminate colonies representing other subspecies (Kauhausen-Keller and Keller 1994). In some cases, beekeepers aim to avoid some subspecies and their hybrids, for example, Africanized honey bees in the American continents (Francoy et al. 2008;Sheppard et al. 1991Sheppard et al. , 1999. There is a range of different methods of honey bee subspecies identification (Bouga et al. 2011), which differ in precision. Usually, methods that are more precise require more effort; therefore, for the purpose of biodiversity conservation and breeding, often rapid and less precise methods were used. Measurement of the cubital index and other traditional methods Samborski et al. 2002;Rostecki et al. 2007) based on a small number of measurements are suitable for the discrimination of a limited number of subspecies. The precision and power of the identification were markedly improved by increasing the number of measured characters (DuPraw 1965a, b). Classical morphometry is based on measurements of 36 characters, which include the size of various body parts (forewings, abdomen and legs), colour and pilosity (Ruttner 1988). However, this method is highly labour intensive because it requires the preparation of several body parts and numerous measurements. Recently, molecular methods based on allozymes (Bouga et al. 2005), mitochondrial DNA (Smith and Brown 1990;Garnery et al. 1993) or microsatellites (Jensen et al. 2005) were used to study the intra-specific variation of honey bees. Among these methods, only microsatellites proved to be suitable for subspecies identification. In general, molecular methods are relatively expensive, require sophisticated equipment and often are more time consuming than morphological methods. Moreover, validation of molecular methods again relies on traditional morphological subspecies identification, and it was demonstrated that morphometric and molecular methods provide similar results (Oleksa and Tofilski 2015).
Efforts have been made to improve and facilitate the identification process based on morphometry. Computerised methods are considered as the most promising (Tofilski 2005;Francoy et al. 2008;Tofilski 2008). In order to reduce the number of dissected body parts, the new methods were only based on forewings and did not require the laborious preparation of sternites, mouthparts and legs. Among other methods of wing measurements, geometric morphometrics proved to be particularly useful. This method was highly effective in the discrimination of three honey bee subspecies used for breeding in Poland (A. m. caucasica , A. m. carnica and A. m. mellifera ) (Tofilski 2008;Gerula et al. 2009). Later, a similar method was used for the discrimination of 24 known honey bee subspecies (Kandemir et al. 2011). Although it was often demonstrated that computerised methods can be useful for the identification of honey bee subspecies, for example, ID-BEES (Batra 1988), FABIS (Rinderer et al. 1990) and ABIS (Francoy et al. 2008), the software is usually not available to the general public. The two notable exceptions are the DrawWing software (Tofilski 2004(Tofilski , 2008 and the ApiClass website (Baylac et al. 2008); however, they cover only a small fraction of honey bee subspecies.
The aim of this study was to develop a fast and easily available method for the identification of a larger number of honey bee subspecies and lineages. The method is based on wing measurements. It requires the identification of 19 points on the forewing image. The method described here was implemented in the IdentiFly software (Tofilski 2017), which is able to identify 19 subspecies and four evolutionary lineages of the honey bee. The software and the data required for the identification are freely available to the general public, including scientists, honey bee breeders and beekeepers. This should make the identification of honey bee subspecies easier and help in the efforts of their conservation.

MATERIALS AND METHODS
We have used images of honey bee worker forewings obtained from the Morphometric Bee Data Bank in Oberursel, Germany. Our dataset was similar to that used in the study of Kandemir et al. (2011) and included 1849 workers from 187 colonies, 25 subspecies and 4 evolutionary lineages (Supplementary Table 1). With very few later additions, this study used the samples, which the analyses in Ruttner's (1988) monography were based on. However, only the part of these samples where wing mounts had been preserved could be included, as the point coordinates used in geometric morphometry cannot be reconstructed from the wing venation angles stored in the original datasets. In order to base our study on a set of well-defined reference samples, we accepted restricted sample numbers in some subspecies, and the wings were re-measured to comply with our current method. We have updated some of the subspecies names used in Ruttner's monography according to the more recent Engel (1999) review (A. m. iberiensis for A. m. iberica , A. m. jemenitica for A. m. yemenitica ) and added A. m. ruttneri (Sheppard et al. 1997). Affiliation of the subspecies to the lineages varies between studies (Ruttner 1988;Garnery et al. 1993;Franck et al. 2001;Kandemir et al. 2011;Meixner et al. 2013). We have used a more recent version of the honey bee taxonomy based in part on genetic studies (Meixner et al. 2013).
Wing pictures were obtained from the Oberursel collection, where they had been taken from mounts on microscopic slides using a camera attached to a stereomicroscope with a resolution of 650 pixels per centimetre (1651 pixels per inch). Each wing was measured by manual identification of 19 characteristic points (later referred to as landmarks, Figure 1). The position of the landmarks was similar to a previous study (Gerula et al. 2009) but was changed slightly for landmarks 7 and 14 in order to agree with that of standard morphometry (Ruttner 1988).
The coordinates of the landmarks were superimpositioned using full Procrustes fit (Dryden and Mardia 1998) in MorphoJ (Klingenberg 2011). Only the aligned coordinates, and not the wing size, were used for statistical analysis in Statistica software v.12 (StatSoft Inc 2011). The analysis was based on average values calculated for each colony. Canonical variate analysis was used for the discrimination of the lineages and subspecies. The correctness of the classification was verified using leave-one-out cross-validation in the PAST 3.11 software (Hammer et al. 2001). All reported classification rates were cross-validated. One of the methods, which can be used for the classification of an unknown individual, is linear discriminant analysis (LDA). In this method, a set of features of the unknown individual is measured, and the measurements are used together with the discriminant functions to calculate the discriminant scores for each of the groups (Fisher 1936). In honey bees, colonies, and not individuals, are usually classified and the mean values of all workers from the colony are used. This allows to achieve higher precision because measurement errors are reduced. The information that is required for classification using this method is the reference configuration and classification functions. Those data are provided as supplementary material (Supplementary Tables 2 and 3) and can be used for classification in both IdentiFly and other software, which might be created in the future. The classification is based on superimposed landmarks. The first step of the classification is Procrustes superposition of the landmarks with the reference configuration. The next step is obtaining discriminant scores for each of the subspecies or lineages. The score is calculated by multiplying the superimposed coordinates with the corresponding discriminant functions and summing them together, including the constant of the discriminant functions. The colony is classified as the subspecies or lineage with the highest classification score. More informative than the classification score is the probability of belonging to a particular class. The probability provided by the IdentiFly software was based on the Mahalanobis distance between an identified colony and the mean of particular classes in the canonical space. For example, if a colony is identified as A. m. mellifera with a probability of 0.01, it means that 1% of colonies in the reference sample of A. m. mellifera is expected to be at a greater distance from the average of this subspecies than the identified colony. In short, a high probability of classification as a particular subspecies means higher similarity of the identified colony to the reference sample of this subspecies. The IdentiFly (Tofilski 2017) software can be downloaded from http://drawwing.org/identifly. It is free for non-commercial use, including research, education and personal use by beekeepers.

Evaluation of the reference database
In the reference data set, the shape of the forewing venation differed markedly between evolutionary lineages (MANOVA F = 34; P < 0.0001). In pairwise comparisons, all lineages differed significantly from each other ( Table I). The largest difference was between lineage C and lineage M ( Table I). The smallest difference was between lineages A and lineage O ( Table I).
Three of the honey bee lineages (A, C and M) were classified correctly without error (100%) (Figure 2). In lineage O, one A. m. anatoliaca colony was incorrectly classified as lineage A. Overall, the cross-validated correct identification rate of lineages was 99.5%.
The shape of the forewing venation also differed significantly among subspecies (MANOVA F = 7; P < 0.0001). In pairwise comparisons, all subspecies differed significantly from each other ( Table II). The largest difference was between A. m. mellifera and A. m. cecropia (Table II). The smallest difference was between A. m. jemenitica and A. m. litorea (Table II). Overall, the cross-validated correct identification rate of subspecies was 88.4%. Nine subspecies were classified without error:  (Table III). Identification rates were high in lineages C and O (91.4 and 97.7%, respectively), where subspecies formed well-separated clusters (Figures 3 and  4). However, they were low (79.7%) in lineage A where most subspecies from central and southern Africa are similar to each other and could not be resolved reliably ( Figure 5).

Program description
To use the IdentiFly software for the identification of lineages and subspecies, forewing images of at least ten workers from one colony need to be obtained as described in the BMaterials and methods^section. Each of the wings has to be measured by indicating 19 landmarks in the appropriate part of the wing. In order to do this, the software should be in the mode BAdd landmarks^. In this mode, clicking the image will result in adding a landmark. Careful positioning of the landmarks (Figure 1) is essential for correct classification. The positions of the landmarks are then saved within the file and can later be verified, corrected and exported to other software. When all wing images from one colony are measured, the colony can be classified. Though it is also possible to classify a single wing from one worker instead of colony samples, this is not recommended because the results are unreliable (Tofilski 2008). In order to properly assign a sample to a subspecies, the classification probabilities, as displayed in a results window, should be taken into account, where a high probability of assignment indicates a high similarity of the classified colony to a particular group of reference samples.
As an example, we present the classification results of one colony collected in Cyprus. First, we have used the lineage classification file and the colony was classified as O lineage (Suppl. Fig. 1). The probability of classification to this group was 0.16, which is much higher than the values for lineages A, C and M: 2.0 × 10 −10 , 1.1 × 10 −8 and 1.5 × 10 −17 , respectively. In this case, it is safe to conclude that the identified colony was from lineage O. In the graph showing CV1 and CV2, the black point representing the colony is within the O ellipse.
Next, we have used the subspecies classification file, and the same colony was classified as A. m. armeniaca (Suppl. Fig. 2). The classification probability (1.5 × 10 −3 ) is much lower than previously. The second highest probability (1.3 × 10 −8 ) is for A. m. syriaca . In this case, we should conclude that the identified colony is similar to A. m. armeniaca but markedly different from the reference sample of this subspecies, and it either is a hybrid with A. m. syriaca or it belongs to a subspecies not covered by the identification model. We know that the latter explanation is true, because by using full standard morphometry, this colony had been identified as A. m. cypria . The colony was correctly classified as belonging to lineage O. However, it could not be classified as A. m. cypria because this subspecies was not included in the classification model due to a small sample size.
As a second example of identification, we have used a colony from Austria. The colony was classified as lineage M (Suppl. Fig. 3) with a probability of 0.42. This number is much higher than the probabilities for other lineages, which were smaller than 10 −11 ; therefore, it is safe to conclude that the classification was correct. When the subspecies classification file was used, the same colony was classified as A. m. mellifera (Suppl. Fig.  4) with a probability of 0.97. The probabilities for all other subspecies were smaller than 10 −12 . In this situation, it can be concluded that the subspecies classification is also correct. However, in lineage       Computer software for identification of honey bee subspecies and evolutionary lineages M, there are two subspecies A. m. mellifera and A. m. iberiensis , with the latter one not included in the classification model due to a small sample size. However, Austria is far outside the natural distribution range of A. m. iberiensis ; thus, the colony is most probably A. m. mellifera , unless it is A. m. iberiensis imported by a beekeeper.

DISCUSSION
The results presented here show that measurements of the forewings of honey bee workers can be successfully used for the identification of evolutionary lineages and subspecies based on the newly developed IdentiFly software. This can be useful for the conservation of endangered subspecies. In a practical context, this may help beekeepers eliminate colonies that differ the most from their preferred native subspecies. The average identification rate of subspecies was close to 90%. This relatively high error margin is mainly due to the A lineage subspecies, where only 79% of colonies were correctly identified. If lineage A subspecies were excluded, the average identification rate would be satisfactorily high (95%). In particular, the identification rates for the commercially used bee subspecies (A. m. mellifera , A. m. carnica and A. m. caucasica ) were satisfactory (Table III). Only the identification of A. m. ligustica is below average (Table III). In comparison to the subspecies, the identification of evolutionary lineages is much more reliable. The discrimination of evolutionary lineages can be particularly useful for the identification of Africanized bees. We did not find a general improvement of subspecies classification results from a twostep procedure, determining the lineage first, and then determining the subspecies within that particular lineage.
The obtained results are consistent with those obtained in earlier studies (Kauhausen-Keller et al. 1997;Kandemir et al. 2011) using similar research material. The earlier studies demonstrated that the identification of the subspecies was possible, but it did not provide enough information for the identification of unknown samples. A lack of this information hindered the research on honey bee biogeography and the conservation of endangered subspecies (Meixner et al. 2013). The  main novelty of this study is to provide data and methods for the identification of honey bee evolutionary lineages and subspecies. The new method is fast and easily available to beekeepers, which should play a key role in the conservation of honey bee subspecies. The alternative method based on standard morphometry (Ruttner 1988) is more time consuming and reference samples of honey bee subspecies are not widely available.
Although the identification of honey bee subspecies based on molecular methods became more available recently, identification based on morphological measurements has an indispensable place. For one, subspecies definitions are historically linked to morphology, and molecular identification needs to be validated against these definitions. Further on, wing morphometry has some advantages over molecular methods. The identification based on wing measurements is relatively easy and fast. In order to obtain the wing images, either a scanner or a digital camera can be used. It is possible to obtain wing images without detaching them from the valuable specimen (Houle et al. 2003), but in most cases, it is easier and faster to dissect them. The equipment required is relatively inexpensive, and nowadays, it is accessible to both scientists and beekeepers.
The number and the position of landmarks used for the wing measurements differed between studies (Meixner et al. 2013). In some studies, the landmark 19 at the distal part of the marginal cell was not included (point 19, Figure 1) (Tofilski 2004). The position of this landmark also differed between studies and it was located either at the front margin of the wing (Kandemir et al. 2011) or at the apex of the marginal cell (Miguel et al. 2011). The position of landmark 14 also varied between studies. Its position was either in the middle of the radius vein (Tofilski 2008;Gerula et al. 2009;Miguel et al. 2011;Kandemir et al. 2011) or at the midpoint between the posterior edge of the radius vein and the anterior edge of the wing (Ruttner 1988). The numbering of the landmarks also varies between studies (Tofilski 2008;Gerula et al. 2009;Miguel et al. 2011;Kandemir et al. 2011). It was suggested that the position and numbering of the landmarks should be standardised in order to facilitate a comparison between studies (Meixner et al. 2013). The landmark numbering scheme presented here may become such a standard because the data can be easily exported to other formats.
The software provides information about the similarity of a colony to different groups. However, it is difficult to make the decision if an investigated colony belongs to a particular subspecies or not. Bees belonging to different subspecies can interbreed and produce hybrids, and there is possibly a range of intermediate phenotypes between the typical phenotypes of Bpure^subspecies. Further, a particular sample may not be a member of any of the subspecies in the data set. The interpretation of the allocation probabilities may give some indications. In an optimal outcome, the probability for one subspecies will be high, and in all remaining subspecies-very low. In this clear case, the focal colony is without doubt classified as that particular subspecies or lineage. However, it is difficult to recommend a fixed threshold value of the classification probability above which an allocation can be considered unequivocal. In another scenario, classification probabilities might be relatively high in two or several subspecies. In this case, the focal colony should be considered as a hybrid between the subspecies with high classification probabilities. Finally, the classification probability can be low in all subspecies. In this case, the classification of the focal colony should be concluded as being doubtful. This can happen in two cases: when the focal colony is very different from all reference samples or when there is an error in the measurements. The first case can occur when the colony belongs to species other than A. mellifera , for example, A. cerana , or the colony belongs to subspecies not included in the classification model, for example, A. m. pomonella . This should not happen very often because the classification presented here covers most of the honey bee subspecies. Low classification probabilities in all subspecies most often occur because there is an error in the measurements, most commonly by swapping the position of two landmarks. It is important to notice that the classification probabilities indicate the similarity of the investigated colony to a reference sample. Unfortunately, the reference sample set of some subspecies is small and may not cover their entire variation range. A colony within the true variation range of a subspecies may thus be assigned only a low similarity to a reference sample.
The current analysis clearly pointed out difficulties with the discrimination of some subspecies. This might partly come from the restriction of the analysis of wing venation (as for example the confusion of A. m. ligustica with A. m. carnica , where colour is a distinctive discriminative body character), but partly also points to deficiencies in the current state of knowledge, as in the African bees, where morphological variation and distribution patterns are notoriously defying distinctive subspecies classifications (Hepburn and Radloff 1998).
A crucial component of the identification tool is the selection of reference samples on which the classification is based. We have used a part of the samples, which were used by Ruttner (1988) in his monography on the taxonomy of honey bees. This first comprehensive description of honey bee subspecies is still widely accepted and still provides a valid framework, which stood the test of time. It thus can be regarded as a lucky circumstance that the original wing mounts of that crucial first morphometric description of the subspecies were still available. Most samples dated back to the times of less intense colony transportation, and quite a number originated from Brother Adam collections (Brother Adam 1983). To this set, we added A. m. ruttneri , but not the more recently described A. m. pomonella (Sheppard and Meixner 2003) and A. m. simensis (Meixner et al. 2011).
Though well documented and preserved in the Oberursel collection and data bank, and covering almost all known subspecies, this reference data set also faces obvious limitations. The most apparent limitation is a low sample size, which compelled us to exclude six subspecies from the set. The estimation of the shape based on geometric morphometrics requires relatively large sample sizes (Cardini and Elton 2007). It is recommended that, in multivariate analysis, the sample size of each group should markedly exceed the number of variables (Arnold 1983). In the case of honey bee wings described by 19 landmarks, the sample size should be markedly larger than 38 colonies, because each landmark consists of two coordinates. This criterion is not met in the case of subspecies. In order to obtain more precise results, the sample size for most of the subspecies should be increased markedly.
A second limitation comes from the unsystematic coverage of the subspecies' geographic range. Within populations believed to belong to one subspecies, some variation occurs, as in black bees (A. m. mellifera ) from Western and Eastern Europe. This problem also calls for an enlargement of the reference data set, but such a crucial step should be considered with utmost care. Currently, there is no general consensus about subspecies reference samples, and the definition of such a set should be based on a more general experts' agreement. A future reference sample set should be large and cover the whole range of a subspecies. This may be facing some difficulties because the exact borders of the distribution may be unclear; the feral populations can be extinct and geographic variation could be affected strongly by beekeeping (De la Rúa et al. 2009). In addition, some of the subspecies may occur in remote and inaccessible geographic regions, or the present subspecies definitions are under debate. In this situation, the preliminary results about subspecies identification were provided here in hope that they will prove to be useful; however, they should be interpreted with some caution until a more expanded and validated reference sample set is available.

CONCLUSIONS
The IdentiFly computer software presented here was shown to be highly effective in the discrimination of the evolutionary lineages of honey bees. It was also satisfactorily effective in discriminating most honey bee subspecies, though it failed in some, particularly within the African A lineage. The reference sample size of many honey bee subspecies is relatively low and needs to be increased in order to achieve higher precision of identification. The software is freely available and can easily be used by both scientists and beekeepers.