Defensive Behavior and Morphometric Variation in Apis mellifera Colonies From Two Different Agro-Ecological Zones of North-Western Argentina

European lineages of Apis mellifera were first introduced into America for beekeeping purposes. A subsequent introduction and accidental release of A. m. scutellata resulted in hybridization events that gave rise to Africanized populations that rapidly spread throughout the continent. In Argentina, Africanized honey bees (AHBs) have been mostly detected in northern regions of the territory, and represent a valuable genetic resource for the selection of stocks with advantageous characteristics for beekeeping. The objective of the present study was to profile honey bee colonies of wild origin with potential beneficial traits for apiculture using morphological, molecular and behavioral traits. Honey bee colonies chosen for evaluation were located in two different agro-ecological regions in north-western Argentina (Tucumán province): The Chaco Depressed Plain (Leales apiary) and the Piedmont (Famaillá apiary). Each apiary was surveyed three times during the 2017–2018 season (mid-season, wintertime, and early spring) for: brood population, phoretic Varroa level and defensive behavior (run, fly, sting, and hang). At the midpoint of the beekeeping season colonies were also characterized by morphometry (45 variables) and mitochondrial haplotypes (COI–COII intergenic region). Apiaries studied showed similar patterns throughout the beekeeping season, for most of the characteristics monitored. However, significant variation in defensive behavior parameters was found between apiaries at the different times of evaluation. Twelve of 45 morphometric variables also showed significant differences between apiaries. The mitochondrial haplotype analysis revealed a high representation of African A4 and A1 haplotypes (91%) in both apiaries. Haplotype variation was associated with morphometric and behavioral traits. Multivariate analyses [principal component analysis (PCA) and principal coordinate analysis (PCoA)] including morphometric and behavior variables explained 65.3% (PCA) and 48.1% (PCoA) of the variability observed between colonies in the first two components. Several morphometric parameters and “fly” behavior were mainly associated with the separation of the colonies. The results from this study point to a possible association between morphometric and behavioral variation and the adaptation of honey bee colonies to differential agro-ecological conditions. We discuss how the detected variation between apiaries can be used for the selection and preservation of honey bee ecotypes in regional breeding programs.

European lineages of Apis mellifera were first introduced into America for beekeeping purposes. A subsequent introduction and accidental release of A. m. scutellata resulted in hybridization events that gave rise to Africanized populations that rapidly spread throughout the continent. In Argentina, Africanized honey bees (AHBs) have been mostly detected in northern regions of the territory, and represent a valuable genetic resource for the selection of stocks with advantageous characteristics for beekeeping. The objective of the present study was to profile honey bee colonies of wild origin with potential beneficial traits for apiculture using morphological, molecular and behavioral traits. Honey bee colonies chosen for evaluation were located in two different agro-ecological regions in north-western Argentina (Tucumán province): The Chaco Depressed Plain (Leales apiary) and the Piedmont (Famaillá apiary). Each apiary was surveyed three times during the 2017-2018 season (mid-season, wintertime, and early spring) for: brood population, phoretic Varroa level and defensive behavior (run, fly, sting, and hang). At the midpoint of the beekeeping season colonies were also characterized by morphometry (45 variables) and mitochondrial haplotypes (COI-COII intergenic region). Apiaries studied showed similar patterns throughout the beekeeping season, for most of the characteristics monitored. However, significant variation in defensive behavior parameters was found between apiaries at the different times of evaluation. Twelve of 45 morphometric variables also showed significant differences between apiaries. The mitochondrial haplotype analysis revealed a high representation of African A4 and A1 haplotypes (91%) in both apiaries. Haplotype variation was associated with morphometric and behavioral traits. Multivariate analyses [principal component analysis (PCA) and principal coordinate analysis (PCoA)] including morphometric and behavior variables explained 65.3% (PCA) and 48.1% (PCoA) of the variability observed

INTRODUCTION
The western honey bee, Apis mellifera (Hymenoptera: Apidae), plays a crucial role in crop pollination and is considered the most important honey producer worldwide Van Engelsdorp and Meixner, 2010;Guzmán-Novoa et al., 2011;Iqbal et al., 2019). This species has shown great adaptive potential, as it has established in diverse environments (Le Conte and Navajas, 2008;Meixner et al., 2013). In their natural range (across Europe, the Middle East, and Africa), more than 26 morphologically and geographically distinct A. mellifera subspecies have been described (Ruttner, 1988).
Using morphology and genetic analysis, A. mellifera subspecies have been assigned to four main evolutionary branches or lineages (A, C, M, and O) (Ruttner, 1988;Franck et al., 2000;Whitfield et al., 2006). A. mellifera ligustica, A. m. mellifera and A. m. carnica subspecies (assigned to C and M European lineages) were established in the Americas for apicultural practices in the early 19 th century (Bierzychudek, 1979;Salizzi, 2014). In 1956, scientists introduced the African subspecies A. m. scutellata (A lineage) into Brazil for the purpose of improving the genetics of honey bees established in tropical climates. An accidental release of these honey bees of African origin and their uncontrolled dispersion, along with the concomitant hybridization with European genotypes, led to a process of "Africanization" of the then resident A. mellifera (Kerr and Nielsen, 1967;Kent, 1988). Africanized populations spread throughout America, increasing the genetic diversity of local resources or ecotypes (Buco et al., 1986;Sheppard et al., 1991;Guzmán-Novoa et al., 2011).
Africanized honey bee (AHB) colonies have retained some characteristic traits from their African ancestors, such as high defensive behavior, tendency to swarm and abscond, tolerance to the mite Varroa destructor, adaptation to subtropical and tropical climates, and a smaller body size (Breed et al., 2004;Schneider et al., 2004;Francoy et al., 2008;Guzmán-Novoa et al., 2011;Rivera-Marchard et al., 2012). Conversely, European honey bees have been associated with low defensiveness, swarming and absconding, along with high honey production and adaptation to temperate climate (De Grandi-Hoffman et al., 1998;Breed et al., 2004;Hunt, 2007;Medina-Flores et al., 2014).
As a result of hybridization processes, behavioral and genetic variation has been observed in Africanized populations from America. This variation is also attributed to the adaptation of honey bees to different geographic and environmental conditions (Southwick and Moritz, 1987;Guzmán-Novoa and Page, 1999;Breed et al., 2004). Several investigations have characterized Africanized populations established throughout America by means of morphology (Buco et al., 1986;Francoy et al., 2008), behavior (Alaux et al., 2009;Guzmán-Novoa et al., 2011), and genetics (Collet et al., 2006(Collet et al., , 2009Whitfield et al., 2006;Acevedo-Gonzalez et al., 2019), including integrative approaches (Rivera-Marchard et al., 2012). In Argentina, a high genetic variability has been detected in commercial and feral colonies (Agra et al., 2018;Calfee et al., 2020). Furthermore, the latter listed authors, confirmed results previously obtained by Sheppard et al. (1999), regarding the presence of populations derived from different lineages: A. mellifera scutellata (A4 and A1), A. mellifera intermissa, and A. m iberiensis distributed mainly in the northern region of the country. The presence of the A1 haplotype suggests that a second influx of honey bees from Africa to South America occurred, and it is considered as another source of Africanization of the Argentinian honey bee populations (Sheppard et al., 1999;Agra et al., 2018).
The honey bee breeding program of Argentina (MeGA, PROAPI) has focused on the selection, preservation and augmentation of honey bee stocks for beekeeping in different agro-ecological regions of the country (Palacio et al., 2000;Bedescarrabure, 2011). The selection criterion includes hygienic behavior, tolerance to brood diseases, low defensiveness, and high productivity. In this framework, the detection and characterization of wild-origin honey bee colonies with desirable characteristics already adapted to the diverse environmental conditions of Argentina is of fundamental importance to the development of sustainable apiculture at the regional level.
In the present study, we describe colony parameters (strength, brood population and phoretic Varroa levels), in conjunction with morphometric, defensive behavior, and genetic characteristics of two apiaries with wild honey bee colonies located in different agro-ecological regions of north-western Argentina (Tucumán province). Colonies were inspected three times (middle-productive season, wintertime, and early spring) during the 2017-2018 season. We discuss our results by considering the potential environmental factors involved in shaping the differences observed between the apiaries evaluated in this study. The results obtained provide a first screening of useful tools for the selection of honey bees with the desirable characteristics of tolerance to colony diseases and low defensiveness while also adapted to a subtropical climate. In addition, this study provides the basis for the need to preserve honey bee Africanized genetic resources in Argentina given the current threat of unpredictable climate change.

Honey Bee Colonies
Apis mellifera colonies from the Leales and Famaillá apiaries were surveyed in the present study. Each apiary was composed of feral colonies that originated from the natural occupation of empty hives by honey bees from surrounding areas. There are no commercial apiaries within a radius of 5 km, where the empty hives were located.
Of the 11 feral colonies selected for the present survey, seven were from FAM and four from LE (Supplementary Table 1). Colonies were selected according to the following criteria: (1) colony strength, more than seven frames covered by honey bees during beekeeping mid-season (hereafter named mid-season), category = 1 based on the "Beekeeping Manual for Subtropical Environments" (Dini and Bedescarrabure, 2011); and (2) colonies naturally-tolerant to the mite V. destructor as determined by survival for one season without mite control treatment. During the survey, the colonies from both apiaries did not receive any acaricide treatment, and were managed using the same protocol of good beekeeping practices (Dini and Bedescarrabure, 2011;Unger et al., 2013).
The inspections of honey bee colonies were performed at three different times during the 2017-2018 season: mid-season (December 2017), wintertime (July 2018), and early spring (September 2018). During the colony survey the following parameters were monitored: colony strength, honey bee brood population, phoretic Varroa, and defensive behavior. At mid-season, 20 nurse bees were randomly taken from the center of each colony and preserved in ethanol 96% (v/v) for further molecular and morphometric analyses (described below).

Colony Status Measurements
Colony strength was assessed by visual inspection of the top of the hive following the procedure described by Unger et al. (2013) and Figini et al. (2017). Three different categories used were as follows: category 1 (at least seven frames covered by honey bees); category 2 (five to seven frames covered by bees); and category 3 (fewer than five frames covered by bees).
The brood population was estimated for each colony as the total area of combs covered by brood according to De Grandi-Hoffman et al. (1998). Briefly, (1) hives were opened and frames sequentially removed; (2) a panel subdivided into quadrants of equal size was superimposed on each frame and an estimate was made of the area covered by brood; and (3) a total count was made of the number of frames fully occupied by brood.
The evaluation of phoretic Varroa was done through the "jar test" according to De Jong et al. (1982) and Dietemann et al. (2013). Briefly, a sample of about 250-300 bees from at least two brood combs (both sides) were swept into a bottle containing 70% v/v ethanol. The bottle was energetically shaken and ethanol filtered through a special mesh and white cloth to separate Varroa mites from the adult worker bees. Once separated, both groups were counted and the percentage of mites present in the sample was calculated based on the total number of bees (% phoretic Varroa = mites/bees × 100) (De Jong et al., 1982;Dietemann et al., 2013).

Genetic Characterization of the Colonies
DNA was extracted from the thorax of one worker from each colony per apiary following the methods outlined in Sheppard et al. (1991). A total of 11 individuals were analyzed. A partial region of the mitochondrial COI-COII intergenic region was amplified using 25µl PCR reactions using primers and conditions described by Hall and Smith (1991) and Lobo Segura (2000) with some modifications according to Agra et al. (2018). The 25 µl PCR reaction mix consisted of 1 µM of each primer, 0.5 mM of PCR nucleotide mix (Genbiotech, Buenos Aires, Argentina), 1.5 mM MgCl 2 (Inbio Highway, Tandil, Argentina), 19 µl reaction buffer (InbioHighway), 1 U Taq Polymerase (InbioHighway), and 5 µl of DNA template. The PCR amplifications were conducted in a MJ PTC-100 thermal cycler (GMI, Ramsey, MN, United States) with a cycling protocol that consisted of an initial denaturation step of 2 min at 94 • C followed by 40 cycles of 30s at 94 • C, 30s at 55 • C, and 1 min at 72 • C, with a final extension step of 2 min at 72 • C. To obtain restriction fragments a10-µl aliquot of each PCR product was digested with HinfI (Promega, Madison, MN, United States) following manufacturer recommendations. Restriction fragments were separated on 4% (w/v) agarose gels, stained with GelRed, and photographed under UV light following Agra et al. (2018).

Morphometric Measurements
Morphometric data of honey bee workers were obtained from the right hind leg, proboscis, and right fore and hind wings. To perform the measurements, the aforementioned body parts were mounted on glass slides (one individual per slide; 10 individuals per colony and apiary), for a total of 110 preparations. The preparations were then photographed, digitized, and morphological characters measured and analyzed.
Geometric morphometry measurements were performed based on the right fore wing, which included 19 homologous, manually plotted, wing vein landmarks ( Figure 2E). The digitized data of right fore wing preparations including the 19 wing landmarks were analyzed using the TPS package (version 1.46, Rohlf, 2010). Thirty-three variables were considered based on their utility and relevance in previously published studies consisting of 32 partial warp (pw) and the centroid size (Francoy et al., 2008). The nomenclature of each measurement was as follows: the term "pw" followed by the letter identifying the axis (x, y) and a number in series (pwx1, pwy1, pwx2, pwy2 . . . pwx16, pwy16).

Defensive Behavior
Parameters of defensive behavior were measured according to the scoring system developed by Ávalos et al. (2014). A score range from 1 to 4 (1 = the lowest intensity of response; 4 = the highest intensity of response) was assigned to the following behavioral parameters: "run" (tendency of worker bees to run on combs), "fly" (tendency of worker bees to fly off the combs during colony manipulations), "sting" (tendency of worker bees to hit the operator's veil), and "hang" (tendency of worker bees to be grouped). To measure the four parameters listed above, hives were opened each time in the presence of the same observer (by applying a minimum amount of smoke per sting) followed by the direct observation of bee behaviors for a period of 30 s.

Statistical Analysis
The values of colony strength were compared between apiaries (LE and FAM) and time (mid-season, wintertime, and early spring) using the Kruskal-Wallis (K-W) test. Post hoc Dunn's test was applied for multiple comparisons. The data of brood population and phoretic Varroa were analyzed by two-way analysis of variance (ANOVA) (fixed factors = Apiary [LE and FAM] and Time [mid-season, wintertime, and early spring]; Apiary × Time). In the case of significant values for the interaction between factors, one-way ANOVAs and post hoc Tukey's HSD tests were performed. The lost colonies were excluded from the analyses.
As a first exploratory analysis of morphological variation between apiaries, bilateral Student's T-tests were separately carried out for all measured morphometric variables. Bonferroni adjustment for multiple comparisons was applied. Subsequently, all morphometric variables were analyzed by two-way ANOVA [fixed factors = Apiary (LE and FAM) and Haplotypes (A1 and A4); Apiary × Haplotype interaction]. C1 haplotype was dropped from the analysis because it was present in only one (LE) of the two evaluated apiaries. Post hoc Bonferroni test was applied for multiple comparisons. Before analyses, the variables were checked for normality with the Shapiro-Wilks test and for homogeneity of variances by Levene's test.
Defensive behavior variables were analyzed by apiary (LE and FAM) and time (mid-season, wintertime, and early spring). The lost colonies were excluded from the analysis. In addition, the same variables were analyzed by haplotype (A1 and A4) at midseason. C1 haplotype was dropped from the analysis because it was present in only one (LE) of the two evaluated apiaries. All analyses were performed using K-W tests. Post hoc Dunn's test was applied for multiple comparisons.
A two-steps analysis was carried out to detect the most informative morphological and defensive behavior variables. First, a principal component analysis (PCA) was performed using all significant variables from the above mentioned statistical tests at mid-season, since all colonies exhibit the highest bee population at this time of the beekeeping season. For morphometric variables we included 10 samples (individuals) per colony. In the case of variables with only one value registered for the colony (defensive behavior) we considered the same value for all individuals from the same colony. A second multivariate analysis, a principal coordinate analysis (PCoA), was run using the Gower coefficient of similarity (Gower, 1971), considering selected variables from the PCA. The selection of morphometry and defensive behavior variables were performed according to their eigenvector coefficients for the two first components of the PCA (variables with a coefficient > 0.25 were selected). All analyses were performed using InfoStat 2016 statistical software (Di Rienzo et al., 2016) as well as SPSS 28.0 version (IBM Corporation 2010).

RESULTS
Honey bee colonies from LE and FAM apiaries (11 colonies in total: four from LE and seven from Famaillá; Supplementary Table 1) were evaluated from middle-productive season (springsummer 2017) to early spring 2018. All colonies in both apiaries showed the highest strength at mid-season (category 1, criterion established for this survey). During wintertime, LE showed 80%  of the colonies (three colonies) in category 1 and 20% (one colony) in category 2, while FAM showed 42.8% of the colonies (three colonies) in category 1, 28.6% (two colonies) in category 2, and 28.6% of the colonies (two colonies) were lost. In early spring, 50% of the colonies (two colonies) in LE apiary were in category 1 and 50% (two colonies) were identified as lost colonies, while for FAM 14.3% of colonies (one colony) were in category 1 and the rest of the colonies were determined as lost (Figure 3 and Supplementary Table 1). Non-significant differences were observed between apiaries (H = 0.05, P = 0.71 K-W test) neither among times of the season (H = 1.75, P = 0.08 K-W test) for colony strength.
Brood population showed similar patterns throughout the season, with no significant difference in the number of frames fully covered by brood between apiaries (F (1 , 17) = 0.0001; P = 0.99; two-way ANOVA). However, significant differences were observed between different times of the season (F (2 , 17) = 26.24; p < 0.001). Specifically, a significantly higher mean number of frames covered with brood was observed in mid-season (5.14 ± 0.27) compared to early spring (4.45 ± 0.52) and wintertime (2.33 ± 0.29) (post hoc comparison Tukey's HSD test) (Figure 4 and Supplementary Table 1). A non-significant interaction between factors (Apiary and Time of the season) was observed (F (2 , 17) = 0.36; P = 0.70).
The percentage of phoretic Varroa appeared to vary throughout the season in both apiaries, but with no significant differences (F (2 , 17) = 1.01; P = 0.38; two-way ANOVA; Figure 5). The mean percentage of phoretic Varroa tended to be higher in mid-season (mean value for both apiaries: 7.05 ± 1.29%) than in early spring (mean value for both apiaries: 3.95 ± 2.52%), with intermediate values during wintertime (mean value for both apiaries: 4.78% ± 1.38%) (Figure 5). Although the values observed in LE seemed to be higher than those in FAM throughout the season, non-significant differences between apiaries were detected (F (1 , 17) = 1.18; P = 0.20). A non-significant interaction was observed between factors (Apiary and Time of the season) (F (2 , 17) = 0.24; P = 0.78).
The genetic determination (COI-COII mitochondrial haplotypes) performed for all colonies from both apiaries (LE and FAM) showed that they were composed of 92% "A" and 8% "C" haplotypes. Within the "A" lineage, 50% corresponded to the A4 haplotype and the other 50% to the A1 haplotype. C1 was the only mitochondrial haplotype present within the "C" lineage (Supplementary Table 1).
A first exploratory analysis showed that 4 of 33 geometric morphometry variables and one of 12 traditional morphometry variables were significantly different between apiaries (Student's T-test, P > 0.05; Supplementary Table 2). Results obtained by the two-way ANOVA (Apiary; Haplotype; Apiary × Haplotype) showed that nine geometric morphometry variables and two traditional morphometry variables differed between apiaries (Table 1A). Specifically, higher values were observed for ProbL and G5 in LE compared with FAM (Bonferroni test; P < 0.05). The remaining traditional and geometric morphometry variables were not significantly different between apiaries.
Morphometric variables also differed between haplotypes. Ten of 33 geometric morphometry variables and two of 12 traditional morphometry variables showed significant differences between haplotypes (A1 and A4) ( Table 1B). The remaining morphometric variables showed no statistically significant results (P > 0.05). For most of the analyzed variables, the twoway ANOVA showed that the Apiary × Haplotype interaction was not significant, thus differences in morphometric variables between apiaries were not influenced by haplotype. However, interactions (Apiary × Haplotype) resulted significant for L5, pwx8, and pwx12 morphometry variables. In these cases, differences in morphometric variables between apiaries depended on the haplotype detected. Specifically, the interaction analysis for L5 variable showed significant differentiation between the two haplotypes from LE apiary [A4 haplotype (LEA4) and A1 haplotype (LEA1)] while LEA4 and FAMA4 showed no significant difference. Pwx8 and pwx12 variables showed significant differentiation between A4 haplotype from LE (LEA4) vs. FAM (FAMA4), and also between LEA4/LEA1 and LEA4/FAMA1 in the case of pwx8 variable (see more details of comparisons in Supplementary Table 3).
Defensive behavior results showed that both apiaries displayed similar patterns for the four measured variables throughout the season (Figure 6). However, the K-W test showed significant differences between apiaries for "fly" (H = 4.89; P = 0.027). When times of the season were compared, significant differences were observed for "fly" and "sting, " specifically, between wintertime and early spring ("fly" [H = 11.83; P = 0.016] and "sting" [(H = 6.73; P = 0.021)]. In addition, border significant differences were observed for "hang" [H = 10.00; P = 0.053]; K-W test) for the same time comparison. LE showed higher values for "fly" and "hang" and lower values for "sting" compared to FAM (Figure 6). Defensive behavior showed non-significant differences for the four variables evaluated between haplotypes (A1 and A4) (P > 0.05; K-W test).
The PCA for combined morphometric and behavior variables showed that the first three components explained 78.3% of the variability among colonies (PC1 43.8%, PC2 21.5%, and PC3 13%; Supplementary Table 4). The distribution of honey bee colonies in the PCA space was mainly explained by the contribution of ProbL, pwy5, pwy6, and "fly, " and associated with higher values on the PC1, while pwx4, pwx8, and pwy11 were main contributors to lower values on the PC1. In relation to PC2, the main positive contribution was due to the centroid size, pwx12 and G5, while lower values were associated with pwx3, pwy11, and pwx12 (Figure 7 and Supplementary Table 4). The distribution of the colonies from each apiary in the plane (X-Y axes) was independent of mitochondrial haplotype (Figure 7 and Supplementary Table 4). Specifically, LE3 (C1 haplotype)  Figure 1). The pattern observed in the scatterplots (PCA and PCoA) showed that colonies belonging to FAM were partially separated from those of LE.

DISCUSSION
The information obtained in the present study allowed us to explore the relevance of morphometric and behavioral variables as indirect indicators of the potential adaptation of AHB populations to subtropical agro-ecological regions of Argentina. Moreover, this study brings valuable information for the characterization and preservation of Africanized populations and supports the need for honey bee breeding programs established at a regional level. Our study reports similar dynamics between apiaries for colony characteristics, mainly adjusted to the high availability of food and similar environmental nectar influx levels. Brood population remains more stable for LE than FAM throughout the beekeeping season, with the latter showing an abrupt drop at the end of the winter. These results in the wintertime could be due to different nutritional status between apiaries for the two agro-ecological regions and associated with food availability or quality. Previous research has proposed that availability of pollen, nectar reserves and quality of stored pollen are the principal reasons for decreases in brood population. Moreover, other studies have shown that A. mellifera can modulate its reproductive rate according to limiting environmental resources, such as availability of food (Le Conte and Navajas, 2008). An example of this adaptive process is the decline of brood as food reserves are depleted. In accordance with this, our results showed colony losses during the winter and early spring, but the absence of colonies presenting an intermediate population size, in the previous inspection, during the mid-season. This could be explained by the strong tendency of AHBs to swarm (Guzmán-Novoa et al., 2011;Uzunov et al., 2014) in the presence of stress factors such as lack of food availability, invasive insects in the hive (such as, ants, and beetles), among other environmental variables.
The colony evaluations performed in the present study showed the absence of significant differences between apiaries with respect to Varroa levels, which were found to be higher than the expected phoretic Varroa values of the region (5% for mite-treated colonies; UDA Los Sarmientos, 2021). This can be attributed to almost all the remaining colonies being naturally tolerant to the mite (with a relatively high load of phoretic Varroa during mid-season), as expected from their Africanized origin. The natural tolerance of AHBs to high Varroa infestation has been previously described (Schneider et al., 2004;Guzmán-Novoa et al., 2011). In line with our results on apiaries located in subtropical climates, Medina-Flores et al. (2014) reported no differences in the levels of infestation between Africanized and pure European colonies from Mexico.
The analysis of morphometric variables showed significant differences between Famaillá and Leales apiaries. However, only a few traditional morphometry characters were variable. A recent study also demonstrated low levels of differentiation among bee ecotypes using traditional wing length morphometry (Calfee et al., 2020). Present results showed that proboscis length could be used to differentiate apiaries. This variation could be linked to the floral species present in the adjoining areas and to the availability of food, and a reflection of the influence of the differential agro-ecological zone where the populations were located. The environmental characteristics of both regions have been previously described. The Depressed Plain region, location of Leales apiary, is an agricultural region represented by the  ProbL, G5, centroid size, pwx3, pwy4, pwy5, pwy6, pwx7, pwx8, pwy11, and pwx12. Behavioral measurements included: "fly." Mitochondrial haplotypes are: A1, A4, and C1. The colonies were named according to the apiary they belong (L: LE and F: FAM) followed by the number of colony -Haplotype, according to the information described in Supplementary Table 1. Ellipses show colonies from the same apiary. Green ellipse: FAM; yellow ellipse: LE. cultivation of sugar cane, soybean and corn with the presence of native plant species of the Chaco Serrano type (Cruzate et al., 2005;Benedetti et al., 2019), while in the Piedmont region, location of Famaillá apiary, major agricultural activity is based on the cultivation of citrus fruits and sugar cane (Benedetti et al., 2019). The latter region is also characterized by the presence of native plant species belonging to the group of Yungas (Cruzate et al., 2005). Currently there is no published information on the characteristics and use of floral resources by honey bees in these regions, hence a detailed characterization of floral shape, structure and composition of each agro-ecological regions of Argentina would be necessary to test the hypothesis of the rapid adaptation of honey bees to their environment. In this regard, previous studies have predicted a rapid dispersion and adaptation of AHBs to new habitats (Cox, 1994;Rivera-Marchard et al., 2012;Ackerman, 2021). A rapid evolution can generate major changes to characteristics such as morphological or behavioral traits, as occurred with the defensive behavior in the gentle AHBs from Puerto Rico (Avalos et al., 2017). Similar results have been described for honey bee ecotypes from Eastern Europe, for which the proboscis length varies across a climatic cline and as a consequence of a broad hybrid zone of A. mellifera subspecies from this region (Alpatov, 1929;Ruttner, 1952;Meixner et al., 2007). Further phenotypic and genomic analysis is necessary to determine if a rapid adaptation of AHB has occurred in South America, specifically in different agro-ecological regions.
Based on the results from our morphometric analyses, geometric morphometry showed greater sensitivity than traditional morphometry to detect significant differences between apiaries. A positive relationship between geometric morphometry variables and apiary location could indicate the presence of an environmental effect, as well as other factors not directly assessed in this study (such as parental gene effects). In addition, we found that sensitivity decreases when analyzing morphometric differences between haplotypes, as non-significant differences were detected, in accordance with recent results by Porrini et al. (2020). These results could be due to phenotypic similarity between A4 and A1 and the ongoing hybridization process present in the northern region of Argentina, and suggest that geometric morphometry would not be sufficiently sensitive as an indirect marker of mitochondrial haplotype for these populations, as proposed by Kandemir et al. (2011) for other populations.
We detected differences in defensive behavior, specifically for "fly" behavior, between apiaries independent of mitochondrial haplotype. We conclude this to likely be the result of paternal and/or environmental effects. It is worth noting that the apiaries being compared in this study had different drone congregation areas. As it has previously been described, this can affect the genetic composition/variability of surrounding apiaries (Collet et al., 2009;Galindo-Cardona et al., 2017. The paternal effect on defensive behavior has also been previously described (De Grandi-Hoffman et al., 1998;Guzman-Novoa et al., 2005). The study of drone contribution to the genetic makeup of apiaries in both agro-ecological regions of Tucuman province is important to document local variations in honey bee defensive behavior as this can be utilized to identify methods and tools for the selection of bee stocks adapted to specific regional environments for breeding programs.
Environmental effects on honey bee defensive behavior have been addressed by several authors (Buco et al., 1986;Southwick and Moritz, 1987;Rivera-Marchard et al., 2012;Nouvian et al., 2016). In particular, Rivera-Marchard et al. (2012) described the importance of climate to the level of defensive behavior response, including the effect of food reserves, and found that the more limited the food, the greater the response. Along similar lines of thought, Scofield and Mattila (2015) proposed that quality and abundance of pollen that honey bee larvae feed on can affect adult behavior. Other authors have described how lack of nutrients and body reserves can affect the health of the colony, which can become vulnerable to Nosema ceranae and V. destructor, increasing colony stress (Invernizzi et al., 2011). Moreover, the impact of agricultural activities in the vicinity of apiaries has been explored, including agro-chemicals as a major factor negatively affecting honey bee populations (Brown and Paxton, 2009). In our study, we observed lower values of the defensive behavior variables during wintering compared to the productive season. This result is likely associated with the combined effect of limited food resources and decreased number of individuals inside the colony during wintertime. Further analyses, considering agricultural activities in surrounding areas and the detection of potential honey bee stressors in the colony assessment will be useful to evaluate the environmental impact on the defensiveness of locally-adapted honey bees in agricultural settings.
The analysis using the COI-COII mitochondrial region was consistent with recently published works (Agra et al., 2018;Calfee et al., 2020;Porrini et al., 2020) that described the presence of genetic variability and a preponderant presence of African lineages (A1 and A4) in feral honey bee colonies mainly established in the north of Argentina. However, a saturation of AHBs in the north of Argentina previously described by Sheppard et al. (1991) was not observed in our study. In addition, Agra et al. (2018) found a predominance of European mitochondrial haplotypes in commercial apiaries, and of A4 and A1 haplotypes in colonies of wild origin from the north of Argentina. Our results on the genetic assessment of wildorigin honey bee colonies are congruent with the mentioned previous studies, as A1 and A4 haplotypes were also detected in a high frequency.
Results from this study revealed significant differences in morphometric and defensive behavior variables between haplotypes, supporting their potential utility for the selection of honey bee stocks with low defensiveness. Furthermore, our results showed that the feral colonies tested have retained some morphometric and defensive behavior characteristics of their African origin, according to the traits previously described by several authors (Breed et al., 2004;Schneider et al., 2004;Francoy et al., 2008;Guzmán-Novoa et al., 2011;Rivera-Marchard et al., 2012).
The combined evaluation of morphometric and defensive behavior variables, performed in our study by multivariate analyses (PCoA and PCA), indicates that maternal lineage (mt haplotype) is not a determining factor explaining the variation among colonies. However, the variation is best explained by colony location (apiary). This could be partially explained by the high hybridization rate of AHB populations in Northern Argentina and the strong influence of the paternal lineage on defensive behavior as previously described by Clarke et al. (2002) and Guzman-Novoa et al. (2005). The multivariate analyses also indicate that some morphometric measures and defensive behavior traits used in this study could be of future utility to track differences in colony response to different agroecological regions.

CONCLUSION
This study determined differences in morphological and behavioral characteristics of honey bees from apiaries located in different agro-ecological zones of Northwest Argentina. The environment, type of farming and agricultural products of each area could potentially impact the nutritional status of A. mellifera colonies and affect access to nutrients from crops and native vegetation. In line with this, we conclude that the morphological and behavioral differences observed between apiaries could be associated with the adaptation of honey bees to specific resources available in the different agro-ecological regions examined.
The study also highlighted the possible role played by feral drone matings in determining the diversity of honey bee population in the regions examined. Further studies of the paternal contribution of the nearest drone congregation areas, using nuclear molecular genetic markers, will bring valuable higher resolution information and guidance to breeding programs of locally adapted honey bee ecotypes.
As a means to support local sustainable apiculture and the preservation of AHBs, this work presents initial findings from a comparison between honey bee populations from different agro-ecological zones of Northwest Argentina and tools that can be used to characterize honey bee ecotypes with desirable characteristics, such as, low defensiveness and tolerance to Varroa.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
EB, GR, MP, and SL conceived and designed the experiments. EB performed the statistical analyzes. MA participated in the genetic analysis. AS participated in the statistical analyses. EB, GG, and LM participated in the survey and contributed to live material. CG participated in the morphometric measurements. EB, AS, and SL wrote the manuscript. All the authors read and accepted the manuscript, contributed to the article, and approved the submitted version.

ACKNOWLEDGMENTS
We would like to thank Miguel Maldonado for his assistance on honey bee colonies management. We are indebted to Alberto Galindo Cardona for his useful suggestions during his stay in EEA Leales and valuable recommendations to EB during her INTA Insertion Fellowship. We thank the editor and the two reviewers for the careful reading of the manuscript and their insightful comments and suggestions. We are indebted to Everett Weber for his helpful and practical recommendations for the data analysis.