© INRA/DIB-AGIB/EDP Sciences, 2010

1. INTRODUCTION

Anatolia is an important gene center of honeybee (Apis mellifera L.) subspecies and ecotypes because of its ecological richness and geographical position (Adam, 1983; Smith et al., 1997). The subspecies found in Anatolia are A. m. caucasica Gorbatchev in the northeast (Adam, 1983; Smith et al., 1997; Guler and Kaftanoglu, 1999a), A. m. anatoliaca Maa in the centre (Bodenheimer, 1942; Adam, 1983; Ruttner, 1988a; Gencer and Firatli, 1999; Guler and Kaftanoglu, 1999b; Kandemir et al., 2000), A. m. meda in the east and southeast Anatolia (Bodenheimer, 1942; Ruttner, 1988a), A. m. carnica Pollmann in the Trace (Smith et al., 1997; Guler and Bek, 2002; Kandemir et al., 2000; Palmer et al., 2000), and A. m. syriaca in the southeast (Bodenheimer, 1942). In addition, some ecotypes have been described at the intersection of the main geographic regions, such as Mugla in the Aegean region and Borçka-Camili in the northeast Anatolia region (Ruttner, 1988a; Guler, 2001) (Fig. 1). Furthermore, the genetic richness of Anatolia has been supported and expressed by recent studies (Bodur et al., 2007; Cesar et al., 2009).

Figure 1 The native honeybee subspecies (A. m. anatoliaca, A. m. caucasica, A. m. carnica and A. m. meda) and the genotype of West anatolia (Aegean) were determined by Ruttner (1988). Samples of this study were collected from which are indicated by on asterisk (*) (Adapted from Smith et al., 1997).

Migratory beekeeping has been widespread in Turkey for three decades. Furthermore, commercial queen bees have been widely available in the past 20 to 25 years (Guler and Demir, 2005). Migratory beekeeping has been prohibited in northeast Anatolia for the past three decades to protect A. m. caucasica However, there has been no other legislation with respect to which subspecies or ecotype is suitable to what region. Queens and colonies reared from different honeybee subspecies (A. m. caucasica, A. m. anatoliaca) have been sold without control for regions with native geographic subspecies. The same is true for other regions and bee subspecies. It is known that native honeybee subspecies and ecotypes might have lost their characteristics because of hybridisation caused by migratory beekeeping, commercial queen bee usage, and uncontrolled mating (Rinderer, 1986; Ruttner, 1988b; Moritz, 1991; Kauhausen-Keller et al., 1997; Lodesani and Costa, 2003; Moritz, 2004). The mixing of the native subspecies has occurred because beekeepers prefer different races for their higher productivity. It is thought that the level of mixing caused by commercial queen-bee usage is particularly high because of the haplodiploid genetic structure of honeybees (Rinderer, 1986; Poklukar and Kezic, 1994); one queen bee can lay about 4 to 5 thousands unfertilised eggs in each season, which all develop into drones (male sexuals). In addition, because they are haploid, one drone can produce 10 million of spermatozoa that are genetically identical to each other. However, there is no information on the effect of imported commercial queens on the other region’s honeybee subspecies, although the value of protecting native races is well recognised (Ruttner, 1988a; Rinderer et al., 1993; Kauhausen-Keller and Keller, 1994). On the other hand, some geographic subspecies have preserved their genetic identifies for a long time in different geographic regions (Whitfield et al., 2006; Delaney et al., 2009). Twenty-four distinct taxonomic groups (subspecies) can be discerned by morphometric methods (Ruttner, 1988a). Therefore, we need additional identification coefficients of many other morphological characteristics for the identification and discrimination of unknown or hybrid worker bee samples.

The aims of this study were: (1) to determine the standard morphometric classification coefficients, canonical discriminant function coefficients, and constant coefficients of the two Turkish native honeybee subspecies (A. m. anatoliaca and A. m. caucasica) most commonly used for commercial queen rearing; and (2) to apply a model for determining the effect of imported commercial queen bees on the morphological characteristics of native honeybees of the Turkeli area of the Sinop province, as well as predicting the origin of any unknown worker bee samples.

2. MATERIALS AND METHODS

2.1. Honeybee material

We studied the morphological characteristics of worker honeybees belonging to two types of beekeeping enterprises: those rearing queens using the original geographic honey bee subspecies and those using commercially reared queens.

2.1.1. Origin of geographic subspecies reference worker bee samples

Samples were obtained in each case from a commercial queen-rearing enterprise. The Caucasian subspecies was evaluated as two populations: Caucasian Ardahan-Posof (CRA-P) and Caucasian Artvin-Camili (CRA-C). Although these are both Caucasian, the honeybee populations of CRA-P and CRA-C are significantly different from each other in morphological, physiological, and behavioural characteristics because of ecological differences between the geographic regions (Guler, 2001; Dodologlu and Genc, 2002; Guler and Alpay, 2005). The Caucasian queen bee subspecies has been reared by the Development Foundation of Turkey (TKV) for more than 25 years in Ankara province and by other enterprises in Artvin (Borçka-Camili) and Ardahan-Posof (CRA-P). Queens of the Anatolia subspecies have been reared by enterprises in the central Anatolia region. In this study, Caucasian samples were collected from Ardahan-Posof (CRA-P, 41° N 42° E, 35 samples) and Artvin-Camili (CRA-C, 41° N 41° E, 35 samples); Anatolian samples (AR-B, 29 samples) were collected from Ankara (Beypazari-Kazan 39° N 32° E) (Fig. 1). Five samples were used from each geographic subspecies to validate the methods.

2.1.2. Worker bee samples from the area importing commercial queen bees

Turkeli, in the Sinop province (41° N 34° E), was selected as the area where commercial queen bees have been used. Turkeli is an isolated area in terms of geography and transport links (Fig. 1). It is located outside migratory beekeeping regions and routes. Adam (1983) described the honeybees of the Sinop (Turkeli and Dikmen) area in the Black Sea region of Turkey and stated that the best honeybees were found in this area. Except for the trade in queen bees, there has been no introduction of honeybees into the area. Trade in queens has taken place for the past 10 to 15 years, although there are villages where no commercial queens have been used. The honeybee samples used for this study were collected from apiaries that have used commercial queen bees for many years and others which have never used them. The worker honeybee samples were collected from six villages in Turkeli: Turhan, Duzler, Yesiloba, Catakguney, Akcabuk, and town centre. A total of 30 samples were randomly selected from five colonies from each village. The distance between the villages from which the samples were collected varied from 3 to 45 km.

2.2. Methods

First, samples of worker native geographic subspecies used for commercial queen rearing were investigated to develop the standards, clustering diagrams, descriptive standard morphological classification function coefficients, canonical discriminant function coefficients, and constant coefficients to be used for determining the subspecies of unknown worker bee samples. The morphological characteristics of 84 worker honeybee samples from CRA-P (30 samples), CRA-C (30 samples), and Anatolia (AR-B, 24 samples) were measured. Multivariate discriminant analyses were used for the determination of the morphological characteristics distinguishing the honeybee samples (Dupraw, 1965; Coley and Lohnes, 1971; Moritz, 1991; Le, 2001). Standard clustering areas of the two original geographic honeybee subspecies in the coordination system were determined.

Second, the morphological characteristics of 30 worker honeybee samples collected from six villages in the Turkeli area, 5 CRA-P, 5 CRA-C and 5 AR-B (total of 45 samples) were investigated. Some samples from Turkeli belonged to the native bees of the area, some were from colonies headed by commercial queens, and some were hybrids. For that reason, in the second step it was not known which colonies were pure or hybrids, the race to which they or their mother queen belonged, or whether they were native or test samples used to verify the method. Therefore, all 45 samples from Turkeli (30 samples) and private geographic subspecies queen-rearing enterprises (15 test samples) were investigated and assumed to be unknown samples. Two score functions were calculated for each sample using these coefficients. Then, with the help of these score functions, the region and clustered area of the unknown samples in the coordination system were calculated and the original sources of the unknown samples were predicted (Coley and Lohnes, 1971; Le, 2001).

2.2.1. Measurement of morphologic characteristics

A total of 129 samples from CRA-P, CRA-C, AR-B, and the Turkeli area collected between June and September of 2006 and 2007 were measured. The samples were stored in ethyl alcohol until preparation. Each part of the samples was fixed on the slide using Hoyer (safe from chloral hydrate) liquid (Borror et al., 1992). Each sample consisted of 15 worker bees from individual colonies, so a total of 1935 (15 × 129) worker bees were used for morphological measurements (Ruttner et al., 1978). Forty-one morphological characteristics were measured as listed in Table I (see also Alpatov, 1929; Goetze, 1940; Dupraw, 1965; Ruttner et al., 1978; Moritz, 1992; Kauhausenkeller and Keller, 1994; Guler and Bek, 2002). Measurements were made using a stereomicroscope. A drawing-tube attachment to the microscope was used for measurement of the vein angels. Eleven vein angels of each wing were measured. First, eighteen different conjugation points of each wing were marked on paper using the ocular microscope, then these points were joined by drawing a line so that the angles between the two lines could be measured (Moritz, 1991; Guler and Bek, 2002).

2.2.2. Calculation of score functionsfor unknown worker bee samplesand prediction of their origins

The effectiveness of the model at predicting the origin of the unknown worker bee samples was determined using 45 samples collected from Turkeli area and three commercial geographic subspecies queen bee rearing enterprises. For that purpose, score functions 1 and 2 were calculated. With the aid of these functions, the clustering areas of unknown samples in the coordinate system were found (Coley and Lohnes, 1971; SAS, 1988). First, the samples were numbered UnS₁, UnS₂, UnS₃, …, UnS₄₅ (unknown samples). This numbering was conducted blindly with regard to the sample source. Score functions 1 and 2 were calculated using the equations below to predict the origin of the unknown samples. To calculate these functions, the standard first discriminant function coefficient (α_i) of each property was multiplied by the phenotypic value of each characteristic (X₁, X₂, X₃, …, X_n). This value was then added to the function 1 constant coefficient (α₀) and the score function 1 calculated (Eq. (1)). Score function 2 was calculated in a similar way (Eq. (2)). In the coordination system, score function 1 is horizontal and score function 2 is vertical (Cooley and Lohnes, 1971). For example, score function 1 and 2 of the sample UnS₁ in Table IV were calculated as explained above and found as shown below. (1)(2)When the score function values were calculated for each of the 45 unknown samples (Tab. IV) and put into the coordination system, Figure 3 was obtained.

2.3. Statistical evaluation

Data were evaluated in two steps. First, multivariate discriminant analysis was applied to the data to find the differences between Caucasian (CRA-P and CRA-C) and Anatolian (AR-B) geographic subspecies in their morphological characteristics and determine the descriptive standard classification function coefficients, discriminant function coefficients, and constant coefficients of the morphological characteristics of two Turkish geographic native subspecies (Le, 2001). Second, the model for the prediction of unknown worker samples was developed using the discriminant function coefficients and constant coefficients of morphological characteristics of these two Turkish native honeybee subspecies (Coley and Lohnes, 1971; SAS, 1988).

3. RESULTS

3.1. Origin of geographic subspecies reference worker bee samples

The mean and standard deviations of the morphological characteristics of Caucasian (CRA-P and CRA-C) and Anatolian (AR-B) worker samples belonging to original queen-rearing enterprises are presented in Table I. There were significant differences between honeybee races in most morphological characteristics, with the exception of WTa, WTb, DWM, WS₆ and the vein angles of J₁₀, L₁₃ and O₂₆.

The first discriminant function described 92.9% of the total variation and discriminated 100% of the two Caucasian ecotypes from the Anatolian race (Wilks’ λ = 0.002). The second discriminant function described 7.19% of the total variation and discriminated the two Caucasian ecotypes (λ = 0.146; Fig. 2). The first discriminant function was correlated (P < 0.01) with LH, LM, WM, WT₃, WT₄, BS, LCb, CT₂, CT₃, CT₄, CSc, and the vein angles of D₇, J₁₆, N₂₃, K₁₉, and O₂₆. The second discriminant function was correlated (P < 0.01) with the other 25 morphological characteristics. All 84 samples (30 CRA-P, 30 CRA-C, 24 AR-B) were classified into their correct groups by the analysis. The groups differed in their distance from the centre. AR-B was located further from the centre than CRA-P + CRA-C, which were located close together in the same direction from the centre (Fig. 2).

3.2. Calculation of morphometric multivariate canonical discriminant function coefficient and constant coefficients

The multivariate canonical discriminant functions of each morphological characteristic and the constant of each function are given in Table II. These coefficients are only valid for the Caucasian (CRA-P and CRA-C) and Anatolian (AR-B) honeybee subspecies measured by this study.

Figure 2 Queen rearing enterprises; analysis of 84 worker honeybee samples from subspecies Caucasian (CRA-P, CRA-C) and Anatolian (AR-B). Horizontal axis: canonical function 1; vertical axis: canonical function 2. Each point represents a sample.

3.3. Worker bee samples from the area importing commercial queens

The mean and standard deviations of the 41 morphological characteristics measured from the 30 worker samples taken from apiaries in Turkeli using commercial queens are presented in Table III. Except for CT₃, D₇, and K₁₉, there were significant differences among worker samples in their morphological characteristics. Samples were significantly different in LCa, J₁₆, and O₂₆ at P < 0.05, different in LH, WT₄, B₄, and E₉ at P < 0.01, and different in all the other characteristics at P < 0.001 levels.

3.4. Calculation of score functions(1 and 2) for unknown worker samples and prediction of their origins

Score function values were calculated for each of the 45 unknown samples (Tab. IV) and put into the coordination system, Figure 3 was obtained, and the origin of the 45 unknown worker bee samples was determined.

Two samples from Turhan, one sample from Duzler, five samples from Yesiloba, four samples from the centre, and five samples from the Akcabuk villages were classified as CRA-P. Two samples from Turhan, one sample from Duzler, and one sample from the centre were classified as CRA-C. One sample from Turhan and three samples from Duzler were classified as AR-B. Five samples from Catakguney (UnS₁₆ UnS₁₇, UnS₁₈, UnS₁₉, and UnS₂₀) were classified as an area different from the original races’ clustered area. The five samples of each of the original test races were all classified correctly. In total, 22 samples (48.89%) out of 45 were classified in CRA-P, eight samples (17.78%) in CRA-C, and 10 samples (22.22%) in AR-B (Fig. 3). Samples from some villages were classified as more than one subspecies. All the samples from Yesiloba and Akcabuk were classified as CRA-P. In addition, score function 1 and 2 were calculated in Microsoft Excel for ease and to standardise the method.

4. DISCUSSION

Multivariate discriminant analysis classified correctly 100% of the 84 original worker honeybee samples taken from native CRA-P, CRA-C, and AR-B geographic honeybee subspecies. Classification level, clustered area, and significance level of the discriminant functions of the races showed that the original geographic subspecies used for queen rearing were from different genetic sources. The high discrimination and classification power showed that the method could accurately predict the origin of unknown worker bee samples. The morphological characteristics of the Caucasian and Anatolian subspecies measured in this study were consistent with those found in previous studies (Guler and Kaftanoglu, 1999a; Gencer and Firatli, 1999; Guler and Kaftanoglu, 1999b; Kandemir et al., 2000; Dodologlu and Genc, 2002). There were also significant differences in morphological characteristics between the two Caucasian ecotypes (CRA-P and CRA-C), as also found in previous studies (Guler, 2001).

Except for the original test samples, almost 57% (17 samples) of the 30 unknown Turkeli samples were clustered in CRA-P and 10% (three samples) in CRA-C, whereas only 16.6% was clustered in the native Anatolian AR-B cluster, showing the impact of Caucasian imported commercial queen usage. Indeed, the Caucasian subspecies has been intensively used in commercial queen-rearing enterprises in Turkey for more than 25 years because it is a calm and productive subspecies. In Turkey, 120000–130000 Caucasian queen bees are reared and sold each year. Recently, the Caucasian CRA-C ecotype has been used for queen rearing. Queen rearing using the Anatolian subspecies is low. These patterns are in accordance with the findings of this study. In addition, the samples classified in CRA-P and CRA-C groups were clustered in the same narrow area and mixed with the original test samples. Therefore, this classification structure and the 100% correct classification of the 15 original test samples proved the reliance of the method used in the study.

Figure 3 Predicted clustering areas of 45 unknown worker honeybee samples using score function 1 and 2. Horizontal axis: canonical function 1; vertical axis: canonical function 2. Each number represents a sample.

Twenty-five out of 30 samples (83.33%) from the Turkeli area, together with the other original test samples, were classified into different clustered areas. This was because of the similarity of the samples in the morphological characteristics measured, the most defining factor for the discriminant functions in predicting the sources of unknown samples. Only five samples (16.67%) were not clustered into the two Turkish native subspecies (Fig. 3), which also had different score function 1 and 2 values from the other 25 unknown samples and 15 original test samples (Tab. IV). These five samples (UnS₁₆, UnS₁₇, UnS₁₈, UnS₁₉, and UnS₂₀) came from the same village (Catakguney) in which there is no use of commercial queens. These samples clustered into a narrow area different from the other samples, showing their homogeneity and that they came from the same genetic source. In addition, these five samples had different morphological characteristics from the original test samples including LPr, LF, LT, LM, WM, WT₃, WT₄, DWM, LFW, LCa, LCb, LS₆, WS₆, vein angles E₉, L₁₃ and J₁₀, LHL, BS, MI, and CSc. Native honeybees of the Turkeli area are different from the Caucasian and Anatolian subspecies because they have a short tongue and a small body with small wings and legs. Catakguney is located in a forested, mountainous area of Turkeli and has only one known beekeeper. Traditional beekeeping practices have been used there, and there has been no colony or queen introduction into this village. All this information supports the result of this study that the honeybees of this village are native bees to the Turkeli area.

Most unknown samples representing the Turkeli area showed similar morphological characteristics to the original test samples clustered in the same group. For instance, all samples from Akcabuk and Yesiloba (10 samples) were classified as CRA-P. Akcabuk samples had the highest mean tongue length (6.642), whereas Yesiloba had the highest vein angle A₄ (34.35°). Ten unknown samples were classified in the AR-B group. According to the score function coefficients, these belonged to the Anatolian race. However, there was a different classification structure in Anatolian clustering area, with the samples forming two groups. The first group had UnS₄₁, UnS₄₂, UnS₄₃, UnS₄₄, and UnS₄₅ and the other had UnS₁, UnS₇, UnS₈, UnS₁₀, and UnS₂₄ (Fig. 2). All the first group samples were original AR-B test samples. They were classified in their real group according to the calculated score functions and clustered into the narrow area near the centre. Although there was no overlap, all the second group samples were clustered close to the CRA-P bees. These samples were collected from three different villages (Turhan, Duzler and the centre) in Turkeli. The samples showed the Anatolian race’s morphological characteristics in HL, TI, LPr, LM, WT₃, WT₄, WS₃, WWM, WS₆, LFW, WFW, and CT₃, as well as vein angels A₄, B₄, D₇, E₉, and N₂₃. For example, these five samples had a 6 to 7 scale value for the colour of tergite 3 (CT₃). This tone of colour is specific to the Anatolian race in Turkey, for which it is scored as 7 to 8 (Ruttner, 1988a; Guler and Kaftanoglu, 1999a). However, these five samples also had some morphological characteristics (CT₂, CT₃, CSc, LHL and BS) specific to the Caucasian honeybee subspecies. Colour of scutellum (CSc) and tergite 4 (CT₃) were dark black rather than black. The Caucasian subspecies has a scale value of 0 to 1 for this characteristic (Bilash et al., 1976; Guler and Kaftanoglu, 1999b). Although some samples showed a similarity with the Caucasian subspecies, others showed a similarity with the Anatolian subspecies or the native bees of Turkeli. This variation was attributed to hybridisation. For that reason, the clustered structure in the AR-B area indicated that some queen bees reared from the breeding colonies of A. m. caucasica but might have mated with drones of A. m. anatoliaca. This finding also shows that hybridisation in any area can be identified using this method.

Although the Turkeli area is outside the migratory beekeeping route and there is no honeybee migration within the Turkeli area, a significantly high morphological variation was found. Catakguney and the centre are 40 km apart, and variation between them might be expected because of this distance. However, there were also morphological differences between the honeybees from Turhan and Duzler, which are only 2 km apart. For instance, two samples from Turhan were classified as CRA-P, two as CRA-C and one as AR-B. Another example is that of Duzler village from which one sample was classified as CRA-P, one as CRA-C, and three as AR-B. There were also differences in WTa, LWM, LFW, WFW, BS, N₂₃, TI and LHL between worker samples collected from the same village, and even from the same apiary. This morphological variation was much higher than found in other studies carried out in Turkey (Guler and Kaftanoglu, 1999a,b; Gencer and Firatli, 1999; Kandemir et al., 2000). Such high morphological variation is not expected for honeybees adapted to a similar geographic region (Alpatov, 1929; Kauhausen-Keller and Keller, 1994; Kauhausen-Keller et al., 1997; Ruttner et al., 2000). These extreme morphological differences suggest that commercial queens reared from the two Turkish original geographic honeybee subspecies have been introduced into Turkeli.

As a result, the model developed also proved that the area is subject to significant genetic mixing. The functions and constant coefficients predicted the origin of the queen bees in 25 of 30 samples of the investigated area, thereby validating the method developed. Morphometric mixing of any original geographic honeybee subspecies and the origin of any unknown worker bee sample can be predicted correctly using the method developed in this study. The validity of the method depends on defining the morphological characteristics of original honeybee subspecies or lines using multiple worker bee samples and correctly determining standard morphometric canonical discriminant function coefficients and constant coefficients.

Acknowledgments

I would like to thank Dr. Yuksel Bek, Ondokuzmayis University Faculty of Medicine and Department of Biostatistics for his statistical advice, and Dr. William Hughes, University of Leeds Institute of Integrative and Comparative Biology for his critical editing of the manuscript. I also thank to Dr DR Smith for giving permission to use the figure adapted from her manuscript.

References

All Tables

All Figures

	Figure 1 The native honeybee subspecies (A. m. anatoliaca, A. m. caucasica, A. m. carnica and A. m. meda) and the genotype of West anatolia (Aegean) were determined by Ruttner (1988). Samples of this study were collected from which are indicated by on asterisk (*) (Adapted from Smith et al., 1997).
In the text

	Figure 2 Queen rearing enterprises; analysis of 84 worker honeybee samples from subspecies Caucasian (CRA-P, CRA-C) and Anatolian (AR-B). Horizontal axis: canonical function 1; vertical axis: canonical function 2. Each point represents a sample.
In the text

	Figure 3 Predicted clustering areas of 45 unknown worker honeybee samples using score function 1 and 2. Horizontal axis: canonical function 1; vertical axis: canonical function 2. Each number represents a sample.
In the text