Introduction

There are at least 24 subspecies of Apis mellifera grouped into four main evolutionary branches: Line O (Near East), line A (African), line M (Western Mediterranean) and line C (Central Mediterranean and Southeastern Europe)¹. In 1956, with the aim of studying adaptation to tropical climates and improving the production of honey, African queen bees (A. mellifera scutellata) and their hybrids were imported from South Africa to Brazil^2,3. In 1957, the escape of bees in Brazil led to massive crossing with local bees of European origin, which generated populations with hybrid genotypes called Africanized⁴. In Mexico, the Africanized bee was first identified in September 1986, when a swarm was caught and identified near Tapachula, Chiapas, in the border area with Guatemala^5,6. Subsequently, the presence of Africanized bees collected in 1987 in the municipalities of Tenosique and Centla, Tabasco, as well as in the municipalities of Coatzacoalcos Veracruz, Palenque Chiapas and Hopelchen Campeche, was reported⁷. In the Yucatan Peninsula, Africanized wild bees (descendants of A. mellifera scutellata) were reported in 1987⁸. Several studies have shown that the bee Apis mellifera scutellata has qualities that deserve to be highlighted; they are mainly honey producers in temperate climates, tolerate cold environments, are prolific suffering from fewer diseases due to parasitosis than other breeds⁹. Africanized bees tend to collect more pollen and more propolis than European bees¹⁰. Due to their greater ability to adapt to the tropical environment, they are widely distributed in the Americas¹¹.

According to FAOSTAT¹², until 2019, China generated the highest production and export of honey globally, with 120,845 t exported, followed by India (65,351 t exported), Argentina (63,522 t), Ukraine (54,834 t), Brazil (30,039 t), Germany (25,239 t) and Mexico (25,122 t). In Mexico, during the period from January to November 2020, honey exports reached 26,077 t, which meant an increase in demand by 3.66 %, compared to 2019. Sixty-two point eighty-six percent of honey production (tonnes) in 2020 was concentrated in seven states, among which Jalisco (6,059), Yucatán (5,529), Chiapas (5,434), Campeche (5,375), Veracruz (4,645), Oaxaca (4,533) and Puebla (2,450) stood out. Tabasco only reported a production of 405 t, ranking 25th as a producer of honey nationwide, much lower than its neighboring states¹³, despite having botanical resources for the development of beekeeping. These data reflect the magnitude of the problem of beekeeping in this entity, located in the humid Mexican tropics.

Currently, beekeeping in Mexico is practiced with Africanized bees and various subspecies of European bees that were introduced to Mexico, such as Apis mellifera mellifera and Apis mellifera ligustica. In the state of Tabasco, the genotypes of bees that coexist, through successive natural crosses, have generated a genetic pool of Africanized bees of which there are not enough reports that allow assessing their population dynamics and their effects on productivity. A limitation is that they are very difficult to identify due to their external morphometric characteristics⁹, however, there are several molecular techniques^{(14, 15)}; which have been successfully used in various studies related to the structure, genetic diversity and phylogeny, determination of mitotypes and gene flow of bee populations^16-22.

In the present study, a comparative analysis between the FABIS (Fast Africanized Bee Identification System) morphometric technique and the RFLP (Restriction Fragment Length Polymorphism) diagnosis was made, using the restriction fragment length polymorphism of mitochondrial DNA (mtDNA) for the identification of Africanized bees from the state of Tabasco and their relationship with 12 wing, femoral and abdominal morphometric variables.

Material and methods

The state of Tabasco is located in the southeast of the Mexican Republic, in the humid tropics region between 17º 15ʹ 00ʺ - 18º 39ʹ 07ʺ N and 90º 50ʹ 23ʺ - 94º 07ʹ 49ʺ W²³.

Results

As can be seen in Table 1, significant statistical differences between the subregions in 7 of the 12 morphometric variables studied (V1 RFWL, V3 RHWL, V5 PRL, V8 FTW, V9 FTL, V10 FTBL and V12 FSL) were found. The bees of the Chontalpa subregion showed on average the largest values in most of the variables, except in V6 TLHL (P= 0.9298), since the largest value was found in the Central subregion (3.02 mm), with a difference of 0.01 and 0.03 mm with respect to the other subregions. No statistically significant differences between the subregions with respect to this variable were found, neither in V2 RFWW (P=0.1368), V4 NHHW (P=0.2764), V7 FLHL (P=0.0945) and V11 FSW (0.5569). The principal component analysis showed that three components have an eigenvalue greater than or equal to 1.0, which together account for 60.088 % of the variability of the original data. CP1 has a positive correlation with all the morphometric variables studied, while CP2 has a positive correlation with five variables and negative or no correlation with the rest. CP3 has a positive correlation with eight variables, and negative or no correlation with the rest.

The graph of influences (Figure 1) shows the coefficients of each variable of the first two components, with right forewing length (V1 RFWL), right hindwing length (V3 RHWL), femur length of the hindleg (V7 FLHL), fourth tergite width (V8 FTW), fourth tergite length (V9 FTL) and fourth sternite length (V12 FSL) being the variables with the greatest influence on CP1. However, the three discriminant variables explain only 60.088 % of the variability of the original data, which is considered insufficient for the purposes of this analysis.

Right forewing length (V1), Right forewing width (V2), Right hindwing length (V3), Number of hamuli of the hindwing (V4), Proboscis length (V5), Tibia length of the hindleg (V6), Femur length of the hindleg (V7), Fourth tergite width (V8), Fourth tergite length (V9), Fourth tergite band length (V10), Fourth sternite width (V11), and Fourth sternite length (V12).

Figure 1 Two-dimensional information of the influences of morphometric variables on two principal components 

According to the diagnosis by the FABIS method, at the state level (n= 153), 67 colonies (43.79 %) were determined with Africanized morphotype, 77 colonies (50.33 %) with European morphotype and 9 colonies (5.88 %) suspicious (S), showing significant differences (Chi²= 52.86, n= 153, P=0.0001). The suspicious colonies were found in apiaries of the Central subregion (n= 31) 1 colony (3.2 %), Chontalpa subregion (n= 68) 6 colonies (8.82 %) and in the Pantanos subregion (n= 24), 2 colonies (8.33 %).

According to the mtDNA analysis (Figure 2, Table 2), at the state level (n= 153), 86 colonies (56.21 %) were determined with Africanized mitotype and 67 colonies (43.79 %) with European mitotype, without finding statistical significance (Chi²= 2.36, n= 153, P=0.124). Specifically, of the 135 commercial colonies sampled, relative frequencies of mitotype E genes were detected in 52 colonies (38.52 %) and 83 colonies (61.48 %) presented mitotype A, finding significant differences (Chi²= 7.12, n= 135, P=0.0076). This demonstrates the predominance of the Africanized mitotype over the European one in an order of 1.59 A/E in the commercial colonies. As can be seen in Table 2, on average, the Ríos subregion showed the lowest frequency of mitotype A (41.67 %), with a higher degree of Africanization found in the rest of the subregions, up to 75 % in the Pantanos subregion.

Af= Africanized, undigested mitotype, Eur= European (double banding). M= molecular size marker.

Figure 2 Amplification of the mtDNA of Apis mellifera digested with the Bgl II restriction enzyme 

When performing the morphometric analysis of the 12 variables using the four categories described as sources of variation, statistically significant differences are found in 67 % of the variables (8 of 12). In Table 3, it can be seen that there are small numerical differences, however, E-BRC bees present in most of the morphometric variables the largest dimensions, except in V5 PRL, however, there are no significant differences.

The clustering analysis presented in Figure 3 shows that, in the first stage of the procedure, category A-CC is grouped with category E-CC (1.63). Next, a second group forms with category E-BRC (it appears as E-CPC in Figure 3) with a grouping distance of 3.96, with category A-WC (it appears as A-CS in Figure 3) appearing as the furthest group (5.76).

Figure 3 Phenogram obtained by the nearest neighbor grouping method of 12 morphometric variables, measured in honeybees (Apis mellifera L.) in the state of Tabasco 

Discussion

By using subregions as sources of variation, the differences found in the variables of length and width of structures, although significant, are very narrow and it is not possible to achieve clear discrimination of bees between subregions. On average, the bees from the state of Tabasco measure V1 RFWL (9.022 mm), V2 RFWW (3.13 mm), V3 RHWL (6.326 mm), V4 NHHW (they have 21.012 hamuli), V6 TLHL (3.01 mm), lower than Creole bees from a population of Africanized bees Apis mellifera sp. from the Lambayeque region in Peru²⁵. However, Tabasco bees are larger in terms of V7 FLHL (2.56 mm).

Pearson’s moment-product correlations between each pair of variables (n= 153) showed that V1 RFWL has a strong positive and significant correlation with V2 RFWW (0.5072), V3 RHWL (0.8368), V7 FLHL (0.5220), V8 FTW (0.5377), V9 FTL (0.5568) and V12 FSL (0.5502), but with the rest the correlations are weak. The non-significant weak correlation (0.1541) between V3 RHWL and V4 NHHW is striking. V5 PRL only has significant positive weak correlation with V9 FTL (0.1754), its correlation with the rest of the variables is weak and not significant. V7 FLHL showed significant correlation with V8 FTW (0.4130), V9 FTL (0.4785), V10 (0.2910), V11 (0.4719) and V12 (0.4550). In relation to the tergite and sternite, the measures correlate positively and significantly with each other. V8 FTW has strong positive and significant correlation with V9 FTL (0.6392), V10 (0.3441), V11 (0.4528) and V12 (0.5667).

Significant differences were found both between the subregions and between the categories, in the variables V1 RFWL and V7 FLHL, which has special attention since they are used for the determination of Africanization by the FABIS method^22,23. Likewise, these variables together showed significant positive correlations with 11 of the 12 variables studied. Regarding the number of hamuli, in a comparative analysis between Africanized and European bees, it was reported that, statistically, there are no differences between Africanized and European worker bees ⁽²⁶; the results of this work also confirm this condition, as the number of hamuli was not found as a discriminating characteristic between both breeds. In relation to the fact that the morphometric pattern of African bees are smaller, and that European bees are larger^3,27,28, under the conditions of this work, it was found that although the morphometric characteristics are different statistically, the difference on average is very narrow, which agrees with several researchers who affirm that Africanized bees are very difficult to differentiate morphometrically, and that currently there are Africanized bees with more European characteristics or vice versa^29,30. This morphometric similarity could lead to erroneous conclusions when using morphometric methods for A-E determination. This was demonstrated in populations of morphometrically Africanized bees³¹ and agrees with these results by finding, in a general way, 59 (38.56 %) non-coincident determinations between both diagnostic methods (FABIS/mtDNA). These results suggest the possibility of finding that bees of different mitotype and collection site share the same morphometric group, so it is not possible to classify the bees from the state of Tabasco by geographical subregion, which could be explained by the uses and customs of distribution of queen bees among producers from different subregions.

The body size of bees has a strong genetic basis, however, it has been shown that the management of colonies, particularly cell size, influences the body size of bees³². Another study that reinforces the idea that beekeeping practices affect the genetic type of bee populations is the one reported in 2007 by Antonio³³, who found that 67.39 % of the colonies in the Comarca Lagunera region were European, which was attributed to the intensity with which beekeepers have carried out the replacement of queens. More recently, other researchers³⁴, through FABIS, reported that 91.49 % of the bee colonies in Mexicali and 67.65 % of the colonies in Ensenada have Africanized morphotypes, however, in the wild colonies, it was found that 100 % in Mexicali and 50 % in Ensenada have Africanized morphotypes, which coincides with what was found in this work, and that could be explained by the differences in the technification of the sampled apiaries.

The fact that, in the Tabasco territory, the bees with mitotype A and E in the commercial colonies have very similar morphometric characteristics could be explained a) by the natural hybridization with the wild colonies since the entry of the African bee; b) by the introduction of commercial queen bees or European breeding bees; c) influenced by the characteristics of the honeycombs. The origin of mitotype E in the state of Tabasco is diverse, since from 2003 to 2010, the ISPROTAB (Institute of Production Systems of the Humid Tropics of Tabasco) acquired inseminated bees of different European breeds from certified farms, to reproduce them in their queen farms and donate F1 queens to producers, in order to promote the annual replacement of queen bees; however, after its cancellation, some beekeepers acquire on their own queen bees from local or foreign farms, which may or may not be certified by SAGARPA (Secretariat of Agriculture, Livestock, Rural Development, Fisheries and Food), with the consequent health risk. So, together with the presence of wild colonies, they are the main source of mitotype A genes, which have a great capacity to multiply, and eventually the fertilization of virgin queens of mitotype E with Africanized drones could be resulting in the proliferation of Africanized genes, since it has been reported that African genes are dominant¹¹.

The results show that apiaries are composed of colonies with both mitotypes in different proportions, which confirms the coexistence of both breed types of bees in some apiaries, and coincides with what was recently reported in a study using mtDNA conducted in seven areas of Buenos Aires, Argentina, for a total of 430 colonies, finding that colonies derived from African bees coexist with European ones in two of the seven areas, in addition to the fact that European mitotypes continue to be more frequent, compared to the results they obtained in 2005³⁵. On the other hand, Quezada-Euan⁸, through analysis of alloenzymes in 25 managed colonies, reported 95 % of AHB (Africanized Honeybee) haplotypes in the states of Chiapas and Tabasco, and 73 % in Yucatan, finding the lowest levels in the states of Michoacán and Jalisco (56 and 40 % respectively), which is very close to what it was found in this work (56.21 %), and they mention that their findings could be related to the intensity of queen replacement practices. In relation to results in this research, was reported 56.21 % of Africanized mitotypes found in commercial colonies in the state of Tabasco, which is below the range reported in a study of Africanization using mtDNA conducted in five populations from the state of Veracruz Mexico and three reference populations, since 60 to 77 % of Africanized mitotype was detected in populations located between 72 and 1,300 masl²⁰.

The two diagnostic techniques used in this work clearly define the breed of the samples of the BRC and WC groups, in which both techniques fully coincide in the results, which is not surprising for the mtDNA technique, knowing its robustness; but they also give evidence that morphometric analysis is capable of defining the breed of bees when these have clearly differentiable or extreme measures; however, when measurements are in ranges shared between breeds, morphometric analysis loses accuracy and results can be erratic, as evidenced by the coincidences and non-coincidences of both diagnostic methods.

In this sense, it is important to mention that the mtDNA analysis manifests the ancestral origin of bees due to the generational stability of mtDNA, but the morphometric analysis reflects the difficulty of defining by breed, since bee populations are interbreeding and have physical and behavioral characteristics of both breeds. In a research carried out in Colombia, without referring to the genotype of bees, it was reported that, in the beekeeping regions of Tolima and Boyacá, the average values of the forewing were from 8.74 mm to 8.63 mm, which indicates African morphological uniformity or very close to the African ones, due to the process of Africanization of the bees from those regions⁹. Similarly, in Colombia, a pattern of amplification of Dra I sites, 16S RFLP, of 87.5 % of African (Apis mellifera scutellata) and 12.5 % European (carnica/ligustica) colonies, was reported, pointing out that Tolima has predominantly Africanized genes^36,37. Regarding the wing length, it was found that in A-WC, A-CC and E-CC bees, they measure on average 8.74 mm, 8.99 mm and 9.08 mm respectively. However, bees with E-CC mitotype from Tabasco have a right forewing length of 9.08 mm, slightly less than those reported in the highlands ≥9.12 mm¹⁹.

Molecular characterization based on mtDNA has become a widely used technique for the study of the differentiation of subspecies or breeds of honeybee Apis mellifera L. by various authors^37-40. This circular and maternally inherited molecule allows characterizing the queen bee through the workers and thus the entire colony and can be considered a marker of the entire colony. Its study has allowed formulating different hypotheses about the evolution of the subspecies of Apis mellifera and defining five evolutionary lineages: lineage A includes the African subspecies such as intermissa and scutellata, lineage M is formed by the subspecies of Western Europe, including mellifera, lineage C formed by the subspecies of Eastern Europe (ligustica or Italian bee), lineage O that includes the subspecies of Near Eastern and the Y that includes the subspecies Apis mellifera yemenitica from Ethiopia¹⁵. Each of these lineages presents a characteristic composition in the sequence of different regions of mitochondrial DNA, as is the case with the subunit I of the cytochrome oxidase (COI) gene, for which lineage M presents a target of the Hinc II endonuclease that is not in lineage A⁴⁰.

Conclusions and implications

Morphometric analyses were only discriminant in the extreme values, and therefore they classify as suspicious of being Africanized, which could be reflecting the reality of the dynamics of intercrossing of populations. The coexistence of bees of European and Africanized type has been affected by some practices such as the replacement of queens, which have substantially influenced the type of bees present in the colonies. The reproductive behavior of honeybees contributes strongly to the maintenance of the levels of Africanization of the colonies, since it allows the introgression of wild genes to the bee populations managed by beekeepers, which can be seen in the commercial bee colonies identified in this work as (CC). Our results clearly and significantly discriminate the category A-ES from the E-BRC, showing the extreme morphometric values; however, the categories A-CC and E-CC showed marked similarities; this wide morphometric range suggests genetic diversity, which must be studied to determine the lineages to which bees belong in this region of the country. So far, the best method to discriminate genotypes E or A is the molecular method.