Phenotypic and genetic characterization of Africanized Apis mellifera colonies with natural tolerance to Varroa destructor and contrasting defensive behavior

Africanized Apis mellifera colonies with promising characteristics for beekeeping have been detected in northern Argentina (subtropical climate) and are considered of interest for breeding programs. Integral evaluation of this feral material revealed high colony strength and resistance/tolerance to brood diseases. However, these Africanized honeybees (AHB) also showed variable negative behavioral traits for beekeeping, such as defensiveness, tendency to swarm and avoidance behavior. We developed a protocol for the selection of AHB stocks based on defensive behavior and characterized contrasting colonies for this trait using NGS technologies. For this purpose, population and behavioral parameters were surveyed throughout a beekeeping season in nine daughter colonies obtained from a mother colony (A1 mitochondrial haplotype) with valuable characteristics (tolerance to the mite Varroa destructor, high colony strength and low defensiveness). A Defensive Behavior Index was developed and tested in the colonies under study. Mother and two daughter colonies displaying contrasting defensive behavior were analyzed by ddRADseq. High-quality DNA samples were obtained from 16 workers of each colony. Six pooled samples, including two replicates of each of the three colonies, were processed. A total of 12,971 SNPs were detected against the reference genome of A. mellifera, 142 of which showed significant differences between colonies. We detected SNPs in coding regions, lncRNA, miRNA, rRNA, tRNA, among others. From the original data set, we also identified 647 SNPs located in protein-coding regions, 128 of which are related to 21 genes previously associated with defensive behavior, such as dop3 and dopR2, CaMKII and ADAR, obp9 and obp10, and members of the 5-HT family. We discuss the obtained results by considering the influence of polyandry and paternal lineages on the defensive behavior in AHB and provide baseline information to use this innovative molecular approach, ddRADseq, to assist in the selection and evaluation of honey bee stocks showing low defensive behavior for commercial uses.

1 Introduction

During the16th century, the European honey bee, Apis mellifera L. (Hymenoptera: Apidae), was introduced into America for beekeeping purposes (1–3). The subspecies brought to the continent by European settlers were Apis mellifera. mellifera L, A. m. ligustica Spinola, A. m. carnica Pollmann, A. m. caucasica Pollmann, A. m. lamarckii Cockerell, A. m. syriaca Skorikov, A. m. cypria Pollmann, and A. m. intermissa Buttel-Reepen (2, 3). These subspecies have continued to be imported ever since, through the commercialization of fertilized queens and queen cells (4, 5). In turn, the first recorded entry of honey bees of African origin, A. m. scutellata Lepeletier, took place in 1956 (6) in the context of a genetic breeding program. Despite the controlled conditions of the program, an accidental escape of these African bees, favored by their colonization capacity and genetic dominance, dispersed in an uncontrolled manner throughout South America by crossing with honey bee populations of European origin, initiating a hybridization process known as “Africanization’’ (7–10). In Argentina, populations derived from African subspecies (A. m. scutellata and A. m. intermissa) and from the Iberian Peninsula (A. m. iberiensis Engel) have been registered and are mainly distributed in the northern region of the country. Particularly, the presence of the last two subspecies of A. mellifera suggests a second source of A and M lineages, respectively (11–13).

The Africanization process marked a milestone in American beekeeping since Africanized honey bee (AHB) colonies retained many advantageous productive traits of their African ancestors, such as active resistance to brood diseases and natural tolerance to pathogens and parasites (14–16), as well as high genetic variability, which favored the emergence of climate-adapted ecotypes (5, 9, 17, 18). However, they also exhibited disadvantageous characteristics such as high defensiveness and swarming tendency (19–23). In addition, AHB exhibited a lower productivity (weight of honey produced) compared with EHB (24, 25).

Given the threat of global warming and the growing need for food equity, AHB colonies are regarded as promising genetic resources for beekeeping, due to their shorter breeding cycle, compared to European honey bees (26, 27), and natural adaptation to subtropical and tropical climates (10, 28, 29). Nonetheless, as approximately 3% of the human population is allergic to honey bee venom (30, 31), the defensive behavior displayed by workers is a topic of interest not only for beekeeping but also for public health (32). Defensive behavior in honey bee colonies is defined as the individual and collective reactions of worker bees in response to external disturbances to the hive (33). It involves different actions such as persuasion (rapid flight, high-pitched buzzing, chasing, hitting, and biting), alarm pheromone release (34–36), and stinging as a direct attack action on the threat (37). This behavior has been described as a highly heritable trait, with genetic dominance and paternal effect (29, 38, 39). The level of defense response of a honey bee colony is related to its sensitivity to the alarm pheromone, visual stimuli, and propensity to sting, run or fly (40), and it requires a sophisticated mechanism of sensory integration, involving olfactory, visual and mechanosensory signals (33, 38). Due to the interest in defensive behavior, genetic studies have been carried out using molecular markers (40, 41), and expression analyses of candidate genes have allowed the identification of molecular pathways associated with the expression of this behavior (42, 43). Age, genetics and adaptation to the environment have been proposed as the main factors associated with defensive behavior in honey bees (44–47). In addition, Gibson et al. (48) have evidenced the existence of asymmetric allelic expression patterns in hybrid honey bees associated with maternal biases and epigenetic regulation (38, 49, 50), which would act as regulators of aggression expression.

Beekeeping is an important economic activity in Argentina, as this country is the fourth world producer and the third exporter of high-quality honeys (51). This commercial activity is mostly performed by small and medium-sized Argentine producers (52, 53). In support of beekeepers, Argentina has implemented honey bee breeding programs to select and multiply bee stocks of European origin with desirable characteristics of production and behavior (54). Particularly, the MeGA program (PROAPI-INTA) focuses on the selection and conservation of local stocks (mostly established in a temperate climate) with demonstrated tolerance or active resistance to mites and brood diseases, low defensiveness, and high productivity (54–56). Conversely, in northern regions of the country the presence of Africanized feral populations with undesirable characteristics for beekeeping (swarming tendency and high defensive behavior) precludes the breeding selection contributing to the genetic background of commercial honey bee colonies (5, 57). Recently, colonies with desirable traits (docility, resistance to brood disease, and adaptation to a subtropical climate) have been characterized in the northwestern region of the country (58, 59), offering an opportunity to preserve honey bee stocks for the development of sustainable apiculture in a subtropical climate. These hybrid colonies have possibly inherited the best characteristics from their European and African origins, as have previously been described in other regions of the American Continent (3, 60, 61).

Next Generation Sequencing (NGS) techniques enable the detection of large numbers of single nucleotide polymorphisms (SNPs) in a cost- and time-efficient manner (62). In A. mellifera, these technologies have facilitated high-resolution genomic studies at the population level, genetic variability and evolution analyses and genotype-phenotype relationship supporting marker-assisted selection for breeding (63–66). Within this group, the Double Digest Restriction-Site Associated DNA (ddRADseq) technique has been widely utilized as an affordable alternative for the genotyping of individuals in model or non-model species with a reduced representation of the genome of interest, thus lowering costs and analysis time (67). Several investigations have been carried out in the last decade using RADseq techniques in A. mellifera (63, 68), including evolutionary history studies (69) and analyses of population structure and variability (70–73).

In this article, we describe the integral characterization of AHB colonies of the same wild origin, located in northwest Argentina (Tucumán province). These colonies possess demonstrated tolerance to the ectoparasitic mite Varroa destructor Anderson and Trueman (Acari: Varroidae) and variability in defensive behavior (58). Evaluations were performed throughout the 2019–2021 productive beekeeping seasons and honey bee colonies with contrasting defensive behavior were analyzed by ddRADseq. The SNP analysis presented here provides a genetic characterization of these contrasting colonies and a set of markers and genomic regions potentially associated with defensive behavior to be further analyzed in breeding programs that seek to improve beekeeping using innovative NGS tools.

2 Materials and methods

2.1 Biological material

We selected a previously characterized A. mellifera colony (named LE2) (58), hereafter named LE0, as a mother colony to generate the biological material analyzed in the present study. The LE0 colony, established in the Leales Apiary (Santa Rosa de Leales, Tucumán, Argentina), represents a hybrid honey bee stock and possesses an Africanized mitochondrial lineage (A1). It has shown high colony strength in spring (category 1), as well as low defensiveness and the ability to survive medium-high levels of V. destructor without acaricide treatment since 2017 (58, 59).

Daughter queens were generated in September 2019 using the traslarvae method (74–76) and nucleus creation according to protocols established for the region (77, 78; Figure 1). Once natural fertilization of the new queens was confirmed, surviving daughter colonies (LE1-9), showing similar high strength status, and the mother colony LE0 were transferred to a new location to facilitate their management. This new apiary was established in Tafí Viejo (26°44’08.8’’S; 65°17’14.3’’W), Tucumán province (Argentina). In this region, the landscape corresponds to the piedmont type (79) with predominant commercial citrus production, followed by sugarcane cultivation (80).

FIGURE 1

Figure 1 Diagram showing the selection process of hybrid A. mellifera colonies for obtaining the biological material for this study. The black arrow indicates selection of colonies; The grey arrow indicates multiplication. Mt. Haplotype: mitochondrial haplotype (genetic characterization). The colonies were surveyed six times during the beekeeping season: late spring-summer 2019, autumn 2020, winter 2020, early spring 2020, 1st productive peak 2020 and 2nd productive peak 2020. Crossed out colonies indicate lost colonies during the time of the survey.

2.2 Colony status measurements

The daughter colonies were established and after three months to strengthen (minimum adult population = at least three frames covered by bees) the survey started. Colony strength (adult population) was estimated by visual inspection of the top of the hive following specific protocols set up for honey bee stocks adapted to a subtropical climate (northern Argentina) (77, 81). The presence of hybrid honey bee colonies (mixed European and African lineages) in northern Argentina has been previously detected and characterized (5, 12, 57) and constitute baseline information for the evaluation proposed in the present study. Briefly, three categories of colony strength were assigned according to the number of hive frames covered with adult honey bees when the top of the beehive is opened: the category 1 is registered when at least seven frames of the hive are covered by bees; category 2 (five to seven frames); and category 3 (four or fewer frames).

The brood population was evaluated for each colony as the total area of combs covered with brood and the number of brood frames (82, 83). For the observations, the bee hives were opened, and the brood frames were sequentially removed. A panel subdivided into quadrants of equal size was superimposed on each brood frame to estimate the average area covered with brood. The number of brood frames fully covered with brood was also registered.

The level of phoretic Varroa was estimated using the “Ethanol wash method” (84, 85) only once during the beekeeping season: (ES) on winter surviving colonies (Figure 1). This method is based on the collection of 250 to 300 worker bees from brood frames in a jar containing 70% v/v ethanol. The bottle is then capped and shaken vigorously to dislodge mites from the bees, and then separated by filtering. The percentage of mites is calculated from the number of bees contained in the sample (% phoretic Varroa = mites/bees*100) (84, 85). This parameter was registered in the analyzed colonies with the purpose of confirming their tolerance to the ectoparasitic mite, V. destructor. All the colonies analyzed in the present study were maintained without acaricide treatment.

2.3 Defensive behavior

Defensive behavior was assessed by opening each honey bee hive with minimal application of smoke and subsequent direct observation for 30 seconds according to Ávalos et al. (86). Four defensive behavior variables were registered: “run” (tendency of worker bees to run on honey bee combs), “fly” (tendency of worker bees to fly during honey bee hive colony manipulation), “sting” (tendency of worker bees to hit the operator) and “hang” (tendency of worker bees to be grouped over the honey bee combs). A score range from 1 to 4 (1 = the lowest intensity of response; 4 = the highest intensity of response) was assigned to each of the four behavior variables measured. All the observations were performed and registered by the same operator. In addition, at each time point, to eliminate the effect of the presence of the alarm pheromones released by guarding honey bee workers (36, 87, 88), the behavior of half of the colonies was measured intercalary, and after a period of 24 h, in the remaining colonies.

A defensive behavior index (DB index) was developed in the present study based on a simple linear model. This index summarizes the weighted values of each defensive behavior variable. In the formula described below, each behavior score (previously registered using the traditional approach described by 86) was multiplied by a fixed numerical value according to its importance for the honey bee hive manipulation by the beekeeper, as follows: 1) The “sting” behavior was considered the greatest impact on the hive management and was, therefore, assigned the highest fixed value per unit (0.15), followed by “fly” (0.05) and “hang” and “run” (0.025) (59).

The DB index values ranged from 0.25 to 1 and were obtained from the following formula:

$$\text{Index DB}=(0.15\text{ x }“\text{sting}”\text{ score})+(0.05\text{ x }“\text{fly}”\text{ score})+(0.025\text{ x }“\text{hang}”\text{ score})+(0.025\text{ x }“\text{run}”\text{ score}).$$

The honey bee colonies were surveyed for the above-mentioned parameters (except phoretic Varroa) six times during the 2019–2021 productive beekeeping season as follows: at the beginning of the beekeeping season during the late spring-summer (December) 2019 (B); autumn (April) 2020 (A); winter (Jun) 2020 (W); early spring (September) 2020 (ES); first Productive Peak (PP1) in spring/summer (December) 2020; and second Productive Peak (PP2) in summer (February) 2021 (Figure 1). The colony surveys were performed during each time of the beekeeping season at peak activity hours of the honey bees (between 10 and 12 am) during sunny days with favorable weather conditions (temperature >20°C). For each survey, defensive behavior and strength category were measured on the same day, the behavioral evaluation was performed first and then the category of the colony.

2.4 Statistical analysis

Colony strength values were compared among daughter (LE2-LE9) and mother (LE0) colonies in the Tafí Viejo apiary at different times of the beekeeping season (B, A, W, ES, PP1 and PP2) using the Kruskal-Wallis (K-W) test. Brood population data was analyzed by a one-way analysis of variance (ANOVA) (factors = Hive [LE] and Time [B, A, W, ES, PP1 and PP2]).

Defensive behavior variables (run, hang, sting, and fly and the CD index) were compared among colonies at different times of the beekeeping season (B, A, W, ES, PP1 and PP2). All analyses were performed using Kruskal-Wallis (K-W) tests. In addition, a multivariate Principal Coordinate Analysis (PCoA) was performed using Gower’s similarity coefficient (89), based on the defensive behavior variables which resulted statistically significant in the above-mentioned analyses and registered at the ES (early spring) time. This specific moment in the beekeeping season is important to determine the defensive behavior of a colony since it is the season when the colony is leaving the wintering and reinitiating the productive stage. Statistical analyses were performed using InfoStat 2016 (90).

2.5 Selection of daughter colonies for NGS analysis

Taking into account the population and behavior parameters registered during the evaluation, the mother (LE0) and two daughter colonies were selected to perform Double Digest Restriction-Site Associated DNA (ddRADseq) according to the following criteria: 1) Stable colony strength parameters during early spring [ES] and productive peaks [PP1 and PP2]), 2) Stable brood population (N° frames with brood >6; % brood/Frame >60) during ES, PP1 and PP2 and, 3) Contrasting and even defensive behavior (most contrasting average values of CD index throughout ES, PP1, PP2) (Figure 1). We have selected colonies that presented even defensive behavior during ES-PP1-PP2, since these are the moments when the most intense productive management takes place and the number of individuals within the colony and the volume of entered food are at maximum values.

With this strategy (selected daughter colonies originated from the same mother colony) we expect to detect phenotypic and genetic differences possibly attributed to the differential paternal origin, according to the polyandry of the species (91, 92). Twenty newly-emerged worker bees were randomly chosen and extracted from a brood frame of each colony under laboratory conditions. Each worker was individually placed in 1.5 ml tubes, frozen in liquid nitrogen for 1 min and then preserved frozen (-80°C) to be processed for NGS analyses.

2.6 DNA isolation

An individual DNA extraction from each preserved worker (whole body) was performed using a modified CTAB Chloroform/octanol protocol (93). First, grinding was performed with liquid nitrogen to powder. The powder was then transferred to a 1.5 ml tube with 900 μl of extraction buffer (1.4 M NaCl, 20 mM EDTA, 100 mM Tris-HCl and sterile H₂0). Tubes were incubated at 65°C for 30 min, maintaining gentle agitation. Then 450 μl of chloroform:isoamyl alcohol (24:1) was added and mixing by inversion was performed for 10 min. The resulting emulsions were centrifuged at 14000 rpm for 30 min. The aqueous phase of each sample was recovered into a clean tube, 5 μl of RNase (10 mg/ml) was added, and the tube was incubated at 37°C for 60 min. Then 450 μl of chloroform:isoamyl alcohol (24:1) solution was added, and the samples were mixed by inversion again. The samples were centrifuged again at 14000 rpm for 30 min to recover the aqueous phase of each sample. Finally, the DNA was precipitated by addition of 600 μl of cold isopropanol (-20°C) and subsequent mixing by inversion. DNA precipitates were recovered by spin centrifugation, the supernatant was discarded, and the pellet was washed with 500 μl of cold 70% ETOH (-20°C). After drying, the precipitates were resuspended in 50 μl of sterile H₂0 (HPLC quality).

The quality of the obtained DNA samples was assessed by electrophoresis in 0.8% w/v agarose gel with GelRed (Biotium) according to the manufacturer’s specifications. Then, DNA samples were sent to the Genomic Unit at IABIMO-CONICET INTA Castelar, Buenos Aires (Argentina), where DNA concentration was measured using Qubit (ThermoFisher). Then, two pools per colony (LE0, LE2 and LE6) were performed with equal amounts of DNA from eight individuals each, yielding a total of six pooled samples (16 workers from each colony).

2.7 Preparation and sequencing of ddRADseq libraries

Double Digest Restriction-Site Associated DNA (ddRADseq) requires a particular combination of two restriction enzymes in the digestion step (63, 67). Regarding the criteria for selecting the enzymes for the A. mellifera genome digestion, EcoRI and MspI had been previously tested with other enzymes involving the ddRADseq technique (68, 71, 73), whereas MboI had been used in mitochondrial DNA studies of the genus Apis (94) and mapping approaches (95, 96). In the present study, two pairs of restriction enzymes were tested for their usage in the digestion step before ddRADseq library construction: the MspI/EcoRI combination, selected according to previous bibliography (73), and the MboI/EcoRI pair, available from the sequencing service (Genomic Unit at IABIMO-CONICET INTA, Argentina). In silico digestion was carried out using the reference genome of A. mellifera Amel_HAv3.1 (NCBI Genome Assembly) and SimRAD package (97) in R. In addition, this test of restriction enzymes was performed using a set of randomly selected DNA samples (obtained in this study) to confirm the size and quantity of the resulting fragments. Pooled DNA samples were processed, and ddRADseq libraries were constructed following the protocol described by Aguirre et al. (98) using the MboI-EcoRI combination in the digestion step. Paired-end sequencing (2x250) was performed on a Novaseq6000 (SAGA-CIMMYT, Mexico).

2.8 Bioinformatic analysis

First, the quality of the batches was checked visually by using FastQC (99). Sequences were filtered by quality using the “process_radtags” program from the Stacks v2.62 package (100) with default parameters. Adapters and poor-quality sequences were removed using Trimmomatic v0.32 (101) with recommended settings. The clean reads were mapped to the reference genome of A. mellifera from NCBI Genome Assembly (Amel_HAv3.1), using Bowtie2 v2.4.5 (102) under default settings by first indexing the genome and then mapping the ddRADseq short reads. To assemble loci according to the alignment positions provided for each read and SNP calling, we ran the “ref_map.pl” program (Stacks). The “populations” program was run afterward to generate population-level summary statistics; the corresponding raw SNP matrix was obtained under Variant Call Format (VCF).

To predict the functional effect of polymorphisms, the Variant Effect Predictor (VEP) from the Ensembl Metazoa (version release 55) (103) was used. The density of SNPs per chromosome was plotted using CMplot package (v.3.1.3) (104) in R (105, 106).

For pathway mapping and functional enrichment analysis, ShinyGO v0.76.3 software (107) was used with default settings, including: pathway database for gene count option was allowed; statistically significant pathways were selected by FDR and sorted by Fold Enrichment. The g:Profiler ve107_eg54_p17_bf42210 was also used with database updated on 15/09/2022 (108). The in-cis potential targets of the SNPs annotated as parts of lncRNAs were searched. Protein-coding genes 10 k to 100 k upstream and downstream of lncRNAs were identified and their function was further investigated using KOBAS (version 3.0) (109) and the NCBI database.

A dendrogram was inferred using the filtered SNPs among the six ddRADseq libraries from the VCF file, based on pairwise identity-by-state (IBS) values using the SNPRelate package version 1.32.2 (110). The dendrogram was constructed using the maximum likelihood hierarchical clustering analysis implemented in the snpgdsHCluster (SNPRelate) option and plotted with the R package ggplot2 (111). Chi-square analyses (112) were performed to visualize the genetic differences between colonies (LE2, LE6 and LE0) in relation to population parameters, namely heterozygosity (HTZ), reference homozygosity (REF HMZ) and alternative homozygosity (ALT HMZ), using InfoStat 2016 software (90).

3 Results

3.1 Colony status analysis

All colonies surveyed showed the greatest strength at the beginning of the survey (category 1, Figure 2; Supplementary Table 1). In autumn, two daughters (LE6, 7) and the mother (LE0) colonies (33%) remained in category 1, while 45% of the colonies (four daughter colonies LE2, 4, 8 and 9) dropped to category 2, and 22% (LE1 and LE3 changed drastically to category 3 (Figure 2). During winter, 67% of the colonies (six daughter colonies LE1, 2, 4, 6, 7 and 8) were in category 2, 22% (2 colonies, LE3 and LE9) in category 3 and only the mother colony (LE0) was in category 1 (Figure 2). In early spring, two colonies (LE1 and LE3) (22%) were considered lost, one colony (LE8) remained in category 2 (11%), and six colonies (the mother colony and five daughters LE2, 4, 6, 7. 8, and 9) were upgraded to category 1 (67%) (Figure 2). During the first productive peak (PP1), 33% (three daughter colonies LE1, 3 and 8) recorded losses and the remaining six colonies (LE2, 4, 6, 7, 9 and the mother colony LE0) were in category 1 (67%). For the second productive peak (PP2) we registered five colonies (daughter LE2, 6, 7, 9 and the mother (LE0) colonies; 55%) in category 1 and 45% of the colonies (daughter colonies LE1, 3, 4, 8) were registered as lost (Figure 2). No significant differences were found for colony strength between colonies (H= 2.04, p=0.96 K-W test). However, a statistically significant difference was observed between the six moments of the beekeeping season evaluated (B, A, W, ES, PP1 and PP2) (H= 25.39, p= <0.0001 K-W test).

FIGURE 2

Figure 2 Colony strength of A mellifera colonies from Tafí Viejo apiary at six specific times during the beekeeping season: late spring-summer 2019 (B), autumn 2020 (A), winter 2020 (W), early spring 2020 (ES), 1st productive peak 2020 (PP1) and 2nd productive peak 2020 (PP2). In different colors, the mean percentage of colonies in category 1 (at least seven frames covered by worker bees), category 2 (five to seven frames covered by worker bees) or 3 (four or fewer frames covered by worker bees) and percentage of lost colonies are shown.

The brood population of the surveyed colonies showed dissimilar patterns throughout the beekeeping season for both the number of brood frames and the percentage of brood per frame, except in early spring. During this time of the season, we observed convergent dynamics of the different colonies in terms of the number of brood frames (Supplementary Table 1). In addition, significant differences were observed for the number of fully covered frames per brood (F=5.83; p<0.0001) and the percentage of brood per frame (F=2.56; p=0.0217) among the colonies. However, non-significant differences were observed among colonies for times of the year. Regarding phoretic Varroa, we observed a high variability among the obtained values for each of the surviving colonies at the time of sampling (ES), with a mean of 4.55% and a deviation equal to 3.89 (Supplementary Table 1).

3.2 Defensive behavior analysis

The analysis of defensive behavior variables (Kruskal-Wallis test) showed significant differences between the mother and daughter colonies for “run” (H=16.00; p=0.0163), “hang” (H=14.76; p= 0.0056), “sting” (H= 18.28; p= 0.0026) and the DB Index (H= 15.16; p=0.0483). However, no significant differences were observed between colonies for “fly” (H= 9.44; p=0.2115). In addition, the same dynamics of defensive behavior was observed between the mother and daughter colonies among the different evaluation times during the beekeeping season (“run” [H=3.14; p=0.5958], “hang” [H=0.68; p= 0.9627], “sting” [H= 4.30; p= 0.3511], DB Index [H= 6.75; p=0.2242] and “fly” [H= 6.78; p=0.1694]) (Supplementary Table 1). The distribution of evaluated honey bee colonies in the two-dimension space (PCoA) confirmed that the behavior variables exhibiting significant differences (run, sting, hang and DB Index) were useful to explain 68.6% of the observed variability in the first two components (PC1 = 46.2; PC2 = 22.4; Figure 3). LE2 and LE6/LE8 and LE7 were located in opposite positions in the X-axis, while LE4 and LE9 and the mother colony (LE0) were centrally located (Figure 3).

FIGURE 3

Figure 3 Principal coordinate analysis (PCoA) scatterplot for A. mellifera colonies based on all defensive behavior variables using Gower’s similarity coefficient (89).

3.3 Selection of daughter colonies and sample preparation for NGS analysis

Two daughter (LE2 and LE6) and the mother (LE0) colonies were selected for NGS assay according to the results registered for their biological and behavior parameters, as follows: the three colonies maintained a stable colony strength (category 1) and stable brood population (N° frames with brood >6; % brood/Frame >60) throughout the ES, PP1 and PP2 time period. Also, the two daughter colonies (LE2 and LE6) displayed contrasting and stable defensive behavior. Specifically, we observed the lowest average CD index (0.31) in the LE2 daughter colony and the highest average CD index (0.70) in the LE6, while the mother colony (LE0) showed a mean value (0.39) throughout the time period (Supplementary Table 1). The in silico enzymatic digestion of the reference DNA (A. mellifera Genome Assembly, Amel_HAv3.1 NCBI) with MboI/EcoRI combination yielded a higher number of fragments (16,035 fragments) of the desired size (~450 to 500 bp) than MspI/EcoRI combination (7,562 fragments). The number of total cut sites in the reference genome was 239,138 for MspI, 63,098 for EcoRI and 61,422 for MboI, hence demonstrating a higher efficiency of the first enzyme combination for the preparation of NGS samples. In addition, both enzyme combinations tested in the laboratory confirmed the in silico results (results not shown).

3.4 ddRADseq analysis

The analysis of the six ddRADseq libraries (LE0, LE2 and LE6; two replicates each) yielded a total of 16,360,032 raw paired-end reads (2x250). After cleaning the raw data, a mean value of 99.4% of reads was retained (Supplementary Table 2). An average alignment rate of 60.5% was obtained after mapping reads to the reference genome. A total of 12,971 variant sites (SNPs) were found to be distributed in 33,053 assembled loci.

Of the 12,971 variant sites processed, 80.40% (10,492 SNPs) were identified as present in protein-coding regions. From this percentage, 73% had a synonymous effect whereas 26% was predicted to cause changes in the amino acid sequence (missense variants). The remaining 13.27% (1,721 SNPs) was cataloged as without description, 4.68% (607 SNPs) as lncRNA, 0.79% (103 SNPs) as miRNA, 0.11% (14 SNPs) as non-translating CDS, 0.08% (10 SNPs) as pseudogenes and 0.66% (86 SNPs) and 0.01% (1 SNP) as tRNA and rRNA, respectively (Figure 4). In addition, 70% of the SNPs that produced a variant effect corresponded to intronic variants, accounting for more than half of the variants detected, while upstream and downstream gene variants accounted for 9% and 7%, respectively (Figure 4). The remaining percentage values were distributed as follows: non-coding transcript variant (5%), intergenic variant (3%), synonymous variant (2%), 3-prime UTR variant (1%), 5-prime UTR variant (1%) and others (3%).

FIGURE 4

Figure 4 Histogram depicting the number of SNPs associated with biological effects according to VEP (Ensembl Metazoa).

The functional annotation of 2,527 genes and subsequent gene enrichment analysis identified 113 GO terms. These terms were classified into functional categories such as cellular component (CC, 12 GO terms), molecular function (MF, 30 GO terms) and biological process (BP, 70 GO terms), as well as KEGGs (Kyoto Encyclopedia of Genes and Genomes) pathways (1 term) (Supplementary Table 3). Enriched GO terms for CC category included plasma membrane, cell periphery, cell junction, intrinsic component of plasma membrane and synapse. The most significantly enriched GO terms for the MF category were associated with DNA-binding transcription factor, transcription regulator, calcium ion binding and G-protein-coupled receptor activities. GO terms enriched in BP were mainly related to the regulation of biological processes, such as nucleic-acid transcription, biosynthesis, signaling, cell communication and developmental processes (See Supplementary Table 3). The neuroactive ligand-receptor interaction category from the KEGG database was shown to be enriched with 31 genes belonging to this term (Supplementary Figure 1).

A total of 20 pathways yielded an FDR value of less than 0.035 (Figure 5). The most enriched pathways included “Neuroactive ligand-receptor interaction” (19 genes), “ABC transporters” (4 genes), “Notch signaling pathway” (4 genes) and “Retinol metabolism” (3 genes). Among the top 20 enriched pathways, “Metabolic pathways”, “Neuroactive ligand-receptor interaction” and “MAPK signaling” together included most of the genes of the data (Figure 5). The SNPs detected in the 5-HT genes are part of the neuroactive ligand-receptor interaction pathway. The identified variants are mostly intronic polymorphisms, potentially involved in the alternative splicing by interfering with splice site recognition (Supplementary Table 4).

FIGURE 5

Figure 5 Significant molecular pathways selected by FDR and sorted by Fold Enrichment generated (ShinyGo tool).

The graphical representation of the SNP distribution (12,971) along the A. mellifera reference genome demonstrated the presence of SNPs in the 16 linkage groups, which was enhanced in LG1 (Figure 6). The statistical analysis of the SNP matrix revealed 142 SNPs differentially present between samples (LE0, LE2 and LE6). Hierarchical cluster analysis of the matrix displayed in a dendrogram (Figure 7; Supplementary Table 5) showed a clear clustering of the samples, with LE0 being more closely related to LE6, which exhibited a higher level of defensive behavior than LE2.

FIGURE 6

Figure 6 SNP density across the 16 chromosomes (linkage groups) of A. mellifera representing the number of SNPs within 1 Mb window size. The horizontal axis represents the chromosome length in Mb. Different colors correspond to SNP density (color scale on the right side of the image). LG, linkage groups.

FIGURE 7

Figure 7 Dendrogram based on hierarchical cluster analysis (pairwise identity-by-state values) from SNP data of six A. mellifera colonies.

From a total of 12,971 polymorphic sites, 647 SNPs were located in protein-coding regions corresponding to 128 genes with mostly known functions. Eighty-two percent of them were annotated as related to different aspects of behavior, namely development, morphology, immune response and caste division, among others. From this SNP set, 128 SNPs were mapped in 21 genes previously reported to be involved in defensive behavior, including dop3 and dopR2, Camkii and Adar, obp9 and obp10 or members of the 5-HT family (Supplementary Table 6).

A population analysis of this set of 128 SNPs revealed that LE0 had the highest percentage of reference homozygosity (REF HMZ) (52.63%) and the lowest percentage of heterozygosity (HTZ) (45.7%). LE6 showed the highest percentage of heterozygosity (55.08%), as well as alternative homozygosity (ALT HMZ) (2.34%) and the lowest percentage of REF HMZ (42.58%). In contrast, LE2 exhibited the lowest value for ALT HMZ (1.56%), together with LE0 (Figure 8). At the statistical level, significant differences were observed between samples (Chi-square Pearson: χ2 = 25.4, p = 0.0046) but not between colonies (χ2 = 5.57, p = 0.2334).

FIGURE 8

Figure 8 Histogram of the percentage of heterozygosity (HTZ), reference homozygosity (REF HMZ) and alternative homozygosity (ALT HMZ) observed for 128 SNPs associated with defensive behavior in the A. mellifera colonies analyzed.

We further investigated the possible role of the lncRNAs, in which 607 SNPs were detected by predicting the potential targets of lncRNAs by in-cis regulation. A total of 1,232 target loci between 10 and 100 kb up and downstream from 178 lncRNAs were predicted in the in-cis role. Among them, 28 genes and 18 microRNAs were included (Supplementary Table 7). Some of these genes included Nxr-1, NLG-5, Dop3, Dnmt1a, Y-y, Silk Fibroin 1 2 3 and 4, E75, nAChRb1, Pla2, mcdp and apamin, among others. Using KOBAS enrichment module, we identified three significantly enriched pathways (p-value<0.05) in the predicted in-cis targets of the lncRNAs. The pathways were related to metabolism, such as alpha-linolenic acid metabolism, oxidative phosphorylation, and retinol metabolism.

4 Discussion

The characterization of A. mellifera ecotypes adapted to subtropical climate and the understanding of their genotypic and phenotypic diversity are essential for efficient breeding, conservation, and management of honey bee germplasm. In the present study, we developed a protocol for the selection of Africanized honey bee stocks based on desirable phenotypic characteristics for beekeeping (colony strength, brood population and defensive behavior) and analyzed the genetic variability of colonies that displayed contrasting defensive behavior as a first approach to understanding possible molecular pathways associated with this behavior in the analyzed honey bee stock.

The results showed homogeneity in the adult population (strength) among the surveyed colonies throughout the beekeeping season. This finding indicates that the mother colony and its daughters had similar dynamics in terms of the number of adult worker bees. However, we found significant differences among colonies at different times of the beekeeping season. These differences were mostly attributed to the difference in the sensitivity of the colony to changes in the environment and to the natural population dynamics of colonies in subtropical climate (58, 113–115) performed by honey bees from each colony during the time surveyed. In addition, we observed significant differences between colonies for both variables of brood population, the percentage of brood per frame and the number of frames with brood. For the percentage of brood per frame, the greatest differences were observed during autumn and early spring, which could be related to the differential response of each colony to the changing climatic conditions and the internal management of colony resources (113, 114). Likewise, the variation registered in the number of brood frames during early spring and the first productive peak could be linked to a differential response to the activation of the natural food input to the colony (33). The reported colony losses in our study might be related to several factors depending on the seasonality. For instance, during early spring, the loss of Africanized colonies is associated with a more marked tendency to swarm due to the reentry of food into the colony and the subsequent explosive population growth (9, 116). During productive peaks, losses are mainly related to avoidance behavior, also highly frequent in Africanized colonies. This behavior is described as a consequence of environmental stressors such as handling, the presence of harmful particles in food sources and excessive rainfall, among others (117).

The experimental design used in the present study, including the analysis of behavior, population, and genetic variables from mother and daughter A. mellifera colonies, revealed that the variability reported for four of the five defensive behavior parameters are potentially associated with genetic differences of paternal origin. Firstly, the environmental factors and beekeeping practices must have affected the colonies in the same way, as they shared location and management. Secondly, according to previous studies (38, 48, 82, 118), defensive behavior displays characteristics of dominance and paternal bias, mainly related to the variability generated by the mating behavior and reproductive system of A. mellifera (reproduction managed by a single fertile female mated by 12 to up to 30 drones; (119, 120). The paternal influence on this character could explain the statistical differences detected in our study between the LE0 colony and the daughter colonies since they share the same maternal background.

The DB (Defensive Behavior) index developed in this study to summarize the defensiveness displayed by the A. mellifera colonies surveyed was consistent with the unweighted values of each of its components. In addition, this index was useful as a selection criterion to detect the two daughter colonies (LE2 and LE6) which exhibited contrasting defensive behavior for ddRADseq analysis. Further surveys under field conditions are needed to test its usefulness as a potential tool for the selection of A. mellifera colonies by low defensive behavior under a subtropical climate.

The ddRADseq technique utilized in the present study set up a new combination of restriction enzymes (MboI/EcoRI) yielding a higher number of DNA fragments of desired size (~450 to 500 bp) than the other enzyme combination (MspI/EcoRI) tested. EcoRI and MspI have been previously tested in A. mellifera, in conjunction with other enzymes involving ddRADseq studies (68, 71, 73). Therefore, the present study provides a precedent for the use of MboI/EcoRI for ddRADseq analysis of Africanized A. mellifera colonies.

Regarding ddRADseq libraries, even though sequence retention after quality filtering was very high, we obtained both a lower alignment rate against the Amel_HAv3.1 reference genome and fewer SNPs between samples in comparison with previous studies (68, 71, 72). We expected these results because we used hybrid genetic material instead of pure subspecies, as shown in the cited studies. These differences can also be attributed to the sample size and the fact that our samples were closely related to each other, being derived from the same mother colony. Further studies with a higher number of colonies should be performed to explore in greater depth the genotypic differences detected.

We identified a total of 12,971 SNPs in the overall comparison of mapped reads against the reference A. mellifera genome distributed in the 16 linkage groups (chromosomes), showing a wide SNP coverage in all linkage groups. In 142 SNPs, we found significant differences between paired comparisons of the analyzed samples. According to the distance matrix generated based on this SNP set, the least defensive daughter (colony LE2) proved to be the most distant from the mother colony, LE0, which was relatively closer to the colony that exhibited the highest defensiveness (LE6). These results made evident the paternal influence on this character and the need for controlled conditions to minimize the defensive behavior in daughter colonies, as previously reported by Giray et al. (46), Büchler et al. (74), Lenoir et al. (121).

The analysis of mapped SNPs in protein-coding regions, evidenced the presence of this type of polymorphisms in 21 genes putatively related to defensive behavior according to previous studies, as modulators of learning and memory, olfactory receptors, neurotransmitters or components of both the nervous and immune systems. Some of these genes were as follows: dop3 and dopR2, from the dopamine family, associated with learning and memory (122); Camkii and Adar, identified as behavior modulators by Whitfield et al. (123); obp9 and obp10 belonging to the OBP family, involved in olfactory sensitivity (124); and members of the 5-HT family (four present among the genes detected in the present study) related to neurological processes through the regulation of octopamine and serotonin receptors (125–127). Recently, Acevedo-Gonzalez et al. (128) have described a direct relationship between the 5-HT family genes and the expression of defensive behavior in Africanized honey bees from Puerto Rico, specifically for genes identified as 5HT2a (alpha and beta). On the whole, the available information should indicate common molecular pathways underlying this behavior in pure and hybrid honey bee lineages.

Notably, in the present study, from the 128 SNP set associated with the defensive response, the most defensive colony, LE6, showed higher alternative homozygosis and heterozygosis than the least defensive daughter colony, LE2, whereas the mother colony had the highest percentage of REF HMZ and the lowest HTZ values. Taking these results into account, we assume that the observed behavioral variation could be closely related to the genetic differences detected by ddRADseq, which have a high paternal influence, in agreement with previous studies (48, 82, 118). The colonies might have different paternal origins, given the polyandry present in the species. These preliminary results need further analysis by surveying a greater number of colonies and generations to support the association between behavior and genetics.

We also detected SNPs in long non-coding RNAs (lncRNAs). Previous studies have confirmed the role of lncRNAs in many biological processes, such as cell differentiation and development, as well as immune responses and tumorigenesis (129). Further analysis of the possible in-cis target regions of the identified lncRNAs performed in the present study revealed that previously studied genes were involved in nervous processes such as Neurexin 1 (Nrx-1), Neuroglin 5 (NLG-5) (130) and D2-like dopamine receptor (Dop3) (131). Other putative target genes include DNA methyltransferase 1a (Dnmt1a) (132), yellow-y (Y-y) (133), ecdysone-induced protein 75 (E75) (134), nicotinic acetylcholine receptor (nAChRb1) (135) and bee venom components such as phospholipase A2 (PLA2), mast cell degranulating peptide (mcdp) and apamin (136–138). This raises the possibility that the lncRNAs found in this study with differential SNP frequencies between sister colonies with contrasting behavior may be acting as regulators of the defensive response and opens the door to further investigations to elucidate the role of these lncRNAs in defensive behavior.

5 Conclusion

Phenotypic and genetic characteristics of an Africanized A. mellifera stock tolerant to V. destructor and displaying variability to defensive behavior were evaluated here. Our results suggest a strong paternal influence, evidenced by the polyandry exhibited by this species (multiple fathers contributing to the expression of this character), in agreement with previous findings. The tools and protocols developed and used in the present study, such as the DB index, MboI/EcoRI restriction enzyme combination and ddRADseq, as well as the SNP dataset generated, could be used in future honey bee selection projects based on defensive behavior and adaptability to a subtropical climate. These promising tools would add innovation to honey bee breeding programs at the regional level, thus improving the productive capacity and competitiveness of small and medium-sized honey producers. The selection and multiplication of honey bee colonies with promising characteristics for beekeeping can be benefited by the saturation of drone congregation zones with desired and known genetic characteristics and/or the use/implementation of reproduction stations and the movement of germplasm as have been previously described (3 and references therein).

Our results set a precedent for the use of ddRADseq for the identification of molecular markers and polymorphisms potentially associated with defensive behavior, mainly those related to the 5HT gene family and lncRNA. The gained information paves the way for further research using a larger number of contrasting daughter colonies for the trait of interest in order to perform linkage analyses between this behavior and major genes involved in its expression.

Data availability statement

The datasets presented in this study can be found at https://www.ncbi.nlm.nih.gov/, accession number: PRJNA1005491.

Author contributions

EB and SL conceived and designed the study. EB, SL, AS, GR and GG participated in the design of experiments and the establishment of colony selection criteria. EB, GG and GR maintained and selected the A. mellifera colonies, EB performed DNA isolation experiments. NA, CVF and CA directed the bioinformatic analyses. EB and CF participated in the bioinformatic analyses. AP and PV processed the DNA samples for sequencing. EB, CF, and SL wrote the manuscript. All the authors read and accepted the manuscript, contributed to the article, and approved the submitted version.

Funding

This study was financially supported by the Instituto Nacional de Tecnología Agropecuaria (INTA) through the projects PD-E4-I079 to SL and PE-E1-I017 to GR and by FONTAGRO project N°: FTG/RF-16112-RG.

Acknowledgments

We would like to thank Miguel Maldonado and Luis Maldonado (EEA Famaillá, INTA) and Enrique Oviedo (IIACS, INTA) for their assistance in the colony management. We are indebted to Paola Fontana and Valentina DiPaoli (EEA Famaillá, INTA) for their assistance in DNA isolation assays. This study has benefited from discussions with MeGA-PROAPI staff (Honey bee breeding program, Argentina). We thank the three Reviewers and the Editor for their insightful comments and suggestions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author AS declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/finsc.2023.1175760/full#supplementary-material

Supplementary Figure 1 | GO enrichment analysis. Neuroactive ligand-receptor interaction category (KEGG).

Supplementary Table 1 | Strength, population and behavior parameters of the analyzed A. mellifera colonies. “-” the parameters were not assessed due to the colony loss. The phoretic varroa were measured for the sole purpose of confirming the status of tolerance of the honey bee stock.

Supplementary Table 2 | Percentage of reads retained after filtering with the “process_radtags” program from Stacks 2 (100).

Supplementary Table 3 | g_GOSt result table containing information about the enriched terms, term sizes and corresponding P-values from gProfiler (108).

Supplementary Table 4 | Pathway enrichment analysis generated with the ShinyGo tool (107).

Supplementary Table 5 | Genetic distances between ddRADseq samples using filtered SNPs and based on identity-by-state (IBS) values (SNPRelate, R package).

Supplementary Table 6 | Annotation, position and population parameters of 128 SNPs mapped in 21 A. mellifera genes and previously reported as involved in defensive behavior.

Supplementary Table 7 | List of lncRNAs IDs from the total set of SNPs (first sheet). LncRNA putative in-cis targets genes 100 kb upstream and downstream (second sheet).

References

1. Ruttner F. Breeding techniques and selection for breeding of the honeybee. British Isles Bee Breeders' Association (1988) (England).

Google Scholar

2. Crane E. The world history of beekeeping and honey hunting. Taylor & Francis, London: Routledge (1999). pp. 720.

Google Scholar

3. Marcelino J, Braese C, Christmon K, Evans JD, Gilligan T, Giray T, et al. The movement of Western honey bees (Apis mellifera L.) among US States and territories: History, benefits, risks, and mitigation strategies. Front Ecol Evol (2022) 10. doi: 10.3389/fevo.2022.850600

CrossRef Full Text | Google Scholar

4. Palacio A, Bedascarrasbure E. Orígenes de las abejas. Rev Super campo (1998) 4(45):23–4.

Google Scholar

5. Calfee E, Agra MN, Palacio MA, Ramírez SR, Coop G. Selection and hybridization shaped the rapid spread of African honey bee ancestry in the Americas. PloS Genet (2020) 16(10):e1009038. doi: 10.1371/journal.pgen.1009038

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Kerr WE, Barbieri MR, Bueno D. Reproduction in African and Italian bees and their hybrids. In 1 deg Congresso Brasileiro Apicultura. Brazil: DEMA, SC, and ACA. (1970) p. 130–5.

Google Scholar

7. Kerr WE, Nielsen RA. Sex determination in bees (apinae). J Apicult. Res (1967) 6:3–9. doi: 10.1080/00218839.1967.11100154

CrossRef Full Text | Google Scholar

8. Kent RB. The introduction and diffusion of the African honeybee in South America. Assoc Pacif. Coast Geogr (1988) 50:21–43. doi: 10.1353/pcg.1988.0009

CrossRef Full Text | Google Scholar

9. Guzmán-Novoa E, Espinoza-Montaño LG, Correa-Benítez A, Guzmán-Novoa G. Colonization, impact and control of Africanized honey bees in Mexico. Vet Mex. (2011) 42:149–78.

Google Scholar

10. Guzmán-Novoa E, Morfin N, de la Mora A, Macías-Macías JO, Tapia-González JM, Contreras-Escareño F, et al. The process and outcome of the africanization of honey bees in Mexico: lessons and future directions. Front Ecol Evol (2020) 404. doi: 10.3389/fevo.2020.608091

CrossRef Full Text | Google Scholar

11. Padilla F, Puerta F, Flores JM, Bustos M. Bees, apiculture and the new world. Arch Zootec. (1992) 41(extra):563–7.

Google Scholar

12. Sheppard WS, Rinderer TE, Garnery L, Shimanuki H. Analysis of Africanized honey bee mitochondrial DNA reveals further diversity of origin. Genet Mol Biol (1999) 22:73–5. doi: 10.1590/S1415-47571999000100015

CrossRef Full Text | Google Scholar

13. Ilyasov RA, Lee ML, Takahashi JI, Kwon HW, Nikolenko AG. A revision of subspecies structure of western honey bee Apis mellifera. Saudi J Biol Sci (2020) 27(12):3615–21. doi: 10.1016/j.sjbs.2020.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Page RE, Guzmán-Novoa E. The genetic basis of disease resistance. Honey bee pests, predators, and diseases. Medina OH:AI Root Co (1997), 469–92.

Google Scholar

15. Mortensen AN, Schmehl DR, Allsopp M, Bustamante TA, Kimmel CB, Dykes ME, et al. Differences in Varroa destructor infestation rates of two indigenous subspecies of Apis mellifera in the Republic of South Africa. Exp Appl Acarol. (2016) 68:509–15. doi: 10.1007/s10493-015-9999-8

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Nganso BT, Fombong AT, Yusuf AA, Pirk CW, Stuhl C, Torto B. Hygienic and grooming behaviors in African and European honeybees—New damage categories in Varroa destructor. PloS One (2017) 12(6):e0179329. doi: 10.1371/journal.pone.0179329

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Buco SM, Rinderer TE, Sylvester HA, Collins AM, Lancaster VA, Crewe RM. Morphometric differences between South American Africanized and South African (Apis mellifera scutellata) honey bees. Apidologie (1986) 18:217–22. doi: 10.1051/apido:19870301

CrossRef Full Text | Google Scholar

18. Sheppard WS, Rinderer TE, Mazzoli JA, Stelzer JA, Shimanuki H. Gene flow between African-and European-derived honey bee populations in Argentina. Nature (1991) 349:782–4. doi: 10.1038/349782a0

CrossRef Full Text | Google Scholar

19. Taylor Jr.OR. The past and possible future spread of Africanized honeybees in the Americas. Bee World (1977) 58(1):19–30. doi: 10.1080/0005772X.1977.11097632

CrossRef Full Text | Google Scholar

20. Kerr WE, Del Rio SDL, Barrionuevo MD. The southern limits of the distribution of the Africanized honeybee in South America [Apis mellifera adansonii]. Am Bee J (1982) 122:196–8.

Google Scholar

21. Spivak M, Gilliam M. Facultative expression of hygienic behaviour of honey bees in relation to disease resistance. J Apic. Res (1993) 32(3-4):147–57. doi: 10.1080/00218839.1993.11101300

CrossRef Full Text | Google Scholar

22. Lobo Segura JA. Highly polymorphic DNA markers in an Africanized honey bee population in Costa Rica. Genet Mol Biol (2000) 23:317–22. doi: 10.1590/S1415-47572000000200013

CrossRef Full Text | Google Scholar

23. Clarke KE, Rinderer TE, Franck P, Quezada-Euán JG, Oldroyd BP. The Africanization of honeybees (Apis mellifera L.) of the Yucatan: a study of a massive hybridization event across time. Evolution (2002) 56(7):1462–74. doi: 10.1111/j.0014-3820.2002.tb01458.x

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Winston ML. The biology and management of Africanized honey bees. Annu Rev Entomol (1992) 37(1):173–93. doi: 10.1146/annurev.en.37.010192.001133

CrossRef Full Text | Google Scholar

25. Rinderer TE, Collins AM. Foraging behavior and honey production. In: The" African" Honey Bee. CRC Press (2019). p. 235–57.

Google Scholar

26. Caron DM. Africanized honey bees in the Americas. Am Bee J (2001) 142:327–8.

Google Scholar

27. Correa-Benítez A, Guzmán-Novoa E. Zootecnia apícola. Introducción a la Zootecnia. México DF: FMVZ-UNAM (2006) p. 403–33.

Google Scholar

28. Guzmán-Novoa E, Page RE Jr. Selective breeding of honey bees (Hymenoptera: Apidae) in Africanized areas. J Econ Entomol (1999) 92(3):521–5. doi: 10.1093/jee/92.3.521

CrossRef Full Text | Google Scholar

29. Breed MD, Guzmán-Novoa E, Hunt GJ. Defensive behavior of honey bees: organization, genetics, and comparisons with other bees. Annu Rev Entomol (2004) 49(1):271–98. doi: 10.1146/annurev.ento.49.061802.123155

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Ewan PW. Venom allergy. BMJ (1998) 316(7141):1365–8. doi: 10.1136/bmj.316.7141.1365

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Burzyńska M, Piasecka-Kwiatkowska D. A review of honeybee venom allergens and allergenicity. Int J Mol Sci (2021) 22(16):8371. doi: 10.3390/ijms22168371

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Bilo BM, Rueff F, Mosbech H, Bonifazi F, Oude-Elberink JNG, EAACI Interest. Group on Insect Venom Hypersensitivity. Diagnosis of Hymenoptera venom allergy. Allergy (2005) 60(11):1339–49. doi: 10.1111/j.1398-9995.2005.00963.x

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Nouvian M, Reinhard J, Giurfa M. The defensive response of the honeybee Apis mellifera. J Exp Biol (2016) 219(22):3505–17. doi: 10.1242/jeb.143016

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Maschwitz UW. Alarm substances and alarm behaviour in social Hymenoptera. Nature (1964) 204:324–7. doi: 10.1038/204324a0

CrossRef Full Text | Google Scholar

35. Collins AM, Rinderer TE, Tucker KW, Sylvester HA, Lackett JJ. A model of honeybee defensive behaviour. J Apic. Res (1980) 19(4):224–31. doi: 10.1080/00218839.1980.11100029

CrossRef Full Text | Google Scholar

36. Wager BR, Breed MD. Does honey bee sting alarm pheromone give orientation information to defensive bees? Ann Entomol Soc Am (2000) 93(6):1329–32. doi: 10.1603/0013-8746(2000)093[1329:DHBSAP]2.0.CO;2

CrossRef Full Text | Google Scholar

37. Millor J, Pham-Delegue M, Deneubourg JL, Camazine S. Self-organized defensive behavior in honeybees. Proc Natl Acad Sci (1999) 96(22):12611–5. doi: 10.1073/pnas.96.22.12611

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Guzmán-Novoa E, Hunt GJ, Page RE Jr., Uribe-Rubio JL, Prieto-Merlos D, Becerra-Guzmán F. Paternal effects on the defensive behavior of honeybees. J Hered (2005) 96(4):376–80. doi: 10.1093/jhered/esi038

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Ellis JD, Hepburn HR. An ecological digest of the small hive beetle (Aethina tumida), a symbiont in honey bee colonies (Apis mellifera). Insectes Soc (2006) 53:8–19. doi: 10.1007/s00040-005-0851-8

CrossRef Full Text | Google Scholar

40. Harpur BA, Kadri SM, Orsi RO, Whitfield CW, Zayed A. Defense response in Brazilian honey bees (Apis mellifera scutellata× spp.) is underpinned by complex patterns of admixture. Genome Biol Evol (2020) 12(8):1367–77. doi: 10.1093/gbe/evaa128

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Lobo NF, Ton LQ, Hill CA, Emore C, Romero-Severson J, Hunt GJ, et al. Genomic analysis in the sting-2 quantitative trait locus for defensive behavior in the honey bee, Apis mellifera. Genome Res (2003) 13(12):2588–93. doi: 10.1101/gr.1634503

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Giray T, Galindo-Cardona A, Oskay D. Octopamine influences honey bee foraging preference. J Insect Physiol (2007) 53(7):691–8. doi: 10.1016/j.jinsphys.2007.03.016

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ávalos A, Fang M, Pan H, Ramirez Lluch A, Lipka AE, Zhao SD, et al. Genomic regions influenpcing aggressive behavior in honey bees are defined by colony allele frequencies. Proc Natl Acad Sci (2020) 117(29):17135–41. doi: 10.1073/pnas.1922927117

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Robinson GE, Page RE Jr. Genetic determination of guarding and undertaking in honey-bee colonies. Nature (1988) 333(6171):356–8. doi: 10.1038/333356a0

CrossRef Full Text | Google Scholar

45. Kolmes SA, Fergusson-Kolmes LA. Stinging behavior and residual value of worker honey bees (Apis mellifera). J N Y. Entomol Soc (1989) 97:218–31.

Google Scholar

46. Giray T, Guzmán-Novoa E, Aron CW, Zelinsky B, Fahrbach SE, Robinson GE. Genetic variation in worker temporal polyethism and colony defensiveness in the honey bee, Apis mellifera. Behav Ecol (2000) 11(1):44–55. doi: 10.1093/beheco/11.1.44

CrossRef Full Text | Google Scholar

47. Hunt GJ. Flight and fight: a comparative view of the neurophysiology and genetics of honey bee defensive behavior. J Insect Physiol (2007) 53(5):399–410. doi: 10.1016/j.jinsphys.2007.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Gibson JD, Arechavaleta-Velasco ME, Tsuruda JM, Hunt GJ. Biased allele expression and aggression in hybrid honeybees may be influenced by inappropriate nuclear-cytoplasmic signaling. Front Genet (2015) 6:343. doi: 10.3389/fgene.2015.00343

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Harrison JF, Hall HG. African-European honeybee hybrids have low nonintermediate metabolic capacities. Nature (1993) 363(6426):258–60. doi: 10.1038/363258a0

CrossRef Full Text | Google Scholar

50. Oldroyd BP, Allsopp MH, Roth KM, Remnant EJ, Drewell RA, Beekman M. A parent-of-origin effect on honeybee worker ovary size. Proc R Soc B: Biol Sci (2014) 281(1775):20132388. doi: 10.1098/rspb.2013.2388

CrossRef Full Text | Google Scholar

51. Manzoni C. Miel Argentina, un negocio que podría ser más dulce (2021). Available at: http://www.senasa.gob.ar/senasa-comunica/noticias/miel-Argentina-de-alta-calidad-endulza-al-mundo#:~:text=Actualmente%20la%20Argentina%20es%20el,de%20apicultores%20y%20de%20colmenas.

Google Scholar

52. Ardanaz S. Analisis Socio-Econoímico de diferentes Modelos de Escala Productiva de Colmenas en la Region Sudoeste de la Provincia de Buenos Aires. (2016). Bachelor Thesis. FCAyF. UNLP.

Google Scholar

53. Arribillaga M. Manejo de la loque americana en Argentina, desde su introducción a la actualidad. Doctoral dissertation. Argentina: Universidad Nacional de la Plata (2021).

Google Scholar

54. Palacio MA, Figini EE, Ruffinengo SR, Rodriguez EM, del Hoyo ML, Bedascarrasbure EL. Changes in a population of Apis mellifera L. selected for hygienic behaviour and its relation to brood disease tolerance. Apidologie (2000) 31(4):471–8. doi: 10.1051/apido:2000139

CrossRef Full Text | Google Scholar

55. Bedescarrabure E. Consolidando la Apicultura como herramienta de desarrollo. Ediciones INTA (2011) 25–39.

Google Scholar

56. Merke J. Dinámica poblacional de Varroa destructor y Apis mellifera L. como herramienta para la selección de abejas tolerantes. Doctoral dissertation, PhD thesis, . Argentina: Universidad Nacional de Mar del Plata (2016).

Google Scholar

57. Agra MN, Conte CA, Corva PM, Cladera JL, Lanzavecchia SB, Palacio MA. Molecular characterization of Apis mellifera colonies from Argentina: genotypic admixture associated with ecoclimatic regions and apicultural activities. Entomol Exp App. (2018) 166(9):724–38. doi: 10.1111/eea.12719

CrossRef Full Text | Google Scholar

58. Bianchi E, Agra MN, García C, Gennari G, Maldonado L, Rodríguez GA, et al. Defensive behavior and morphometric variation in Apis mellifera colonies from two different agro-ecological zones of north-Western Argentina. Front Ecol Evol (2021) 9:590225. doi: 10.3389/fevo.2021.590225

CrossRef Full Text | Google Scholar

59. Bianchi E. Bases genéticas del comportamiento defensivo de colonias de Apis mellifera adaptadas a clima subtropical. pH Thesis. Facultad de Ciencias Naturales e Instituto “Miguel Lillo”. Argentina: Universidad Nacional de Tucumán (2022).

Google Scholar

60. Carneiro FE, Torres RR, Strapazzon R, Ramírez SA, Guerra JCV Jr., Koling DF, et al. Changes in the reproductive ability of the Varroa destructor (Anderson & Trueman) in Africanized honey bees (Apis mellifera L.) (Hymenoptera: Apidae) colonies in southern Brazil. Neotropical. Entomol (2007) 36:949–52. doi: 10.1590/S1519-566X2007000600018

CrossRef Full Text | Google Scholar

61. Galindo-Cardona A, Acevedo-Gonzalez JP, Rivera-Marchand B, Giray T. Genetic structure of the gentle Africanized honey bee population (gAHB) in Puerto Rico. BioMed Cent Genet (2013) 14:65. doi: 10.1186/1471-2156-14-65

CrossRef Full Text | Google Scholar

62. Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, et al. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PloS One (2011) 6(5):e19379. doi: 10.1371/journal.pone.0019379

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Stolle E, Moritz RF. RESTseq–efficient benchtop population genomics with RESTriction Fragment SEQuencing. PloS One (2013) 8(5):e63960. doi: 10.1371/journal.pone.0063960

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Grozinger CM, Robinson GE. The power and promise of applying genomics to honey bee health. Curr Opin Insect Sci (2015) 10:124–32. doi: 10.1016/j.cois.2015.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Guarna MM, Hoover SE, Huxter E, Higo H, Moon KM, Domanski D, et al. Peptide biomarkers used for the selective breeding of a complex polygenic trait in honey bees. Sci Rep (2017) 7:8381. doi: 10.1038/s41598-017-08464-2

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Saelao P, Simone-Finstrom M, Avalos A, Bilodeau L, Danka R, De Guzman L, et al. Genome-wide patterns of differentiation within and among U.S. commercial honey bee stocks. BMC Genomics (2020) 21:704. doi: 10.1186/s12864-020-07111-x

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Peterson BK, Weber JN, Kay EH, Fisher HS, Hoekstra HE. Double digest RADseq: An inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PloS One (2012) 7(5):e37135. doi: 10.1371/journal.pone.0037135

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Petersen GE, Fennessy PF, Van Stijn TC, Clarke SM, Dodds KG, Dearden PK. Genotyping-by-sequencing of pooled drone DNA for the management of living honeybee (Apis mellifera) queens in commercial beekeeping operations in New Zealand. Apidologie (2020) 51:545–56. doi: 10.1007/s13592-020-00741-w

CrossRef Full Text | Google Scholar

69. Wallberg A, Han F, Wellhagen G, Dahle B, Kawata M, Haddad N, et al. A worldwide survey of genome sequence variation provides insight into the evolutionary history of the honeybee Apis mellifera. Nat Genet (2014) 46(10):1081–8. doi: 10.1038/ng.3077

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Noel G, Chaichoompu K, Leclercq G, Fabre E, Gengler N, Steen K, et al. (2017). Human based genetic tools to refine genetic populations structure of honey bees (Apis mellifera ssp L.) (Hymenoptera: Apoidea: Apidae) colonies at regional scale, in: Conference: EUROPEAN PhD NETWORK "INSECT SCIENCE" VIII Annual Meeting, Naples, Italy, November 2017.

Google Scholar

71. Petersen GEL, Fennessy PF, Van Stijn TC, Clarke SM, Dearden PK. Genotyping-by-sequencing for genetic improvement in honeybees. In Proc. Assoc. Advmt. Anim. Breed Genet (2017) 22:217–20.

Google Scholar

72. Eimanifar A, Kimball RT, Braun EL, Ellis JD. Mitochondrial genome diversity and population structure of two western honey bee subspecies in the Republic of South Africa. Sci Rep (2018) 8(1):1333. doi: 10.1038/s41598-018-19759-3

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Vaudo AD, Fritz ML, López-Uribe MM. Opening the door to the past: accessing phylogenetic, pathogen, and population data from museum curated bees. Insect Sys. Div (2018) 2(5):4. doi: 10.1093/isd/ixy014

CrossRef Full Text | Google Scholar

74. Büchler R, Andonov S, Bienefeld K, Costa C, Hatjina F, Kezic N, et al. Standard methods for rearing and selection of Apis mellifera queens. J Apic. Res (2013) 52(1):1–30. doi: 10.3896/IBRA.1.52.1.07

CrossRef Full Text | Google Scholar

75. Simbaña Chorlango H. Evaluación de tres métodos de reproducción de abejas reinas de la especie (Apis mellífera) en el cantón Pedro Moncayo 2012 (Bachelor's thesis). (2015).

Google Scholar

76. Oré JC, Sotelo A, Martos A, Chura J. Tres tipos de colmenas relacionado a la crianza y el desarrollo biológico de reinas Apis mellifera. Anales Científicos (2020) 81, 1:266–77. doi: 10.21704/ac.v81i1.1636

CrossRef Full Text | Google Scholar

77. Dini CB, Bedascarrasbure EL. Manual de apicultura para ambientes subtropicales: Una propuesta de la Red de Escuelas del Noroeste Argentino (NOA). (2011).

Google Scholar

78. Figini EE, Barreto JA. (2017). Available at: https://inta.gob.ar/documentosdocumentos/apicultura-multiplicacion-de-colmenas#:~:text=La%20reproducci%C3%B3n%20de%20la%20colonia,n%C3%A9ctar%2C%20polen%20y%20equilibrio%20sanitario.

Google Scholar

79. Collantes M, Busnelli J. El Cuaternario de la provincia de Tucumán. Geología Tucumán INSUGEO (2014), 157–69.

Google Scholar

80. Morales CC, Moreno RA, Benedetti PE, Sanchez HA. Campanãa citriícola 2020: Tucumaín cuenta con 50.472 ha cultivadas (2020). Available at: https://inta.gob.ar/noticias/campana-citricola-2020-tucuman-cuenta-con-50472-ha-cultivadas.

Google Scholar

81. Unger N, Poffer D, Frigoli L, Marcó O, Fourquet G. Manual De Prácticas Apícolas Para Producir Miel de Calidad en la Cuenca Del Salado. Argentina: INTA (2013). Available at: https://inta.gob.ar/sites/default/files/script-tmp-inta_-_manual_de_prcticas_apcolas.pdf.

Google Scholar

82. DeGrandi-Hoffman G, Collins A, Martin JH, Schmidt JO, Spangler HG. Nest defense behavior in colonies from crosses between Africanized and European honey bees (Apis mellifera L.)(Hymenoptera: Apidae). J Insect Behav (1998) 11:37–45. doi: 10.1023/A:1020862432087

CrossRef Full Text | Google Scholar

83. Solignac J, Delpiano J, Arza R, Figin E, Spagnuolo C, Poffer D, et al. Evaluación de suplementos proteicos en colonias de Apis mellifera. In: Memoria técnica de la 1ra reunión anual de la INTA. Buenos Aires Argentina: Estación Experimental Agropecuaria General Villegas (2014). Available at: https://inta.gob.ar/sites/default/files/inta_mt2014_solignac_evaluacion_suplementos_proteicos.pdf.

Google Scholar

84. De Jong D, Roma DDA, Goncalves LS. A comparative analysis of shaking solutions for the detection of Varroa jacobsoni on adult honeybees. Apidologie (1982) 13(3):297–306. doi: 10.1051/apido:19820308

CrossRef Full Text | Google Scholar

85. Dietemann V, Nazzi F, Martin SJ, Anderson DL, Locke B, Delaplane KS, et al. Standard methods for varroa research. J Apic. Res (2013) 52(1):1–54. doi: 10.3896/IBRA.1.52.1.09

CrossRef Full Text | Google Scholar

86. Ávalos A, Rodríguez-Cruz Y, Giray T. Individual responsiveness to shock and colony-level aggression in honey bees: evidence for a genetic component. Behav Ecol Sociobiol. (2014) 68:761–71. doi: 10.1007/s00265-014-1

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Bortolotti L, Costa C. Chemical Communication in the Honey Bee Society. In: Mucignat-Caretta C, editor. Neurobiology of Chemical Communication. Boca Raton (FL: CRC Press/Taylor & Francis (2014).

Google Scholar

88. López-Incera A, Nouvian M, Ried K, Müller T, Briegel HJ. Honeybee communication during collective defence is shaped by predation. BMC Biol (2021) 19:106. doi: 10.1186/s12915-021-01028-x

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Gower JC. A general coefficient of similarity and some of its properties. Biometrics (1971) 27:857–71. doi: 10.2307/2528823

CrossRef Full Text | Google Scholar

90. Di Rienzo JA, Casanoves F, Balzarini MG, Gonzalez L, Tablada M, Robledo CW. InfoStat versión 2016. Grupo InfoStat, FCA. Córdoba, Argentina: Universidad Nacional de Córdoba (2016).

Google Scholar

91. Kraus FB, Neumann P, Van Praagh J, Moritz RFA. Sperm limitation and the evolution of extreme polyandry in honeybees (Apis mellifera L.). Behav Ecol Sociobiol (2004) 55:494–501. doi: 10.1007/s00265-003-0706-0

CrossRef Full Text | Google Scholar

92. Oldroyd BP, Fewell JH. Large fitness benefits from polyandry in the honey bee, Apis mellifera. TREE (2008) 23(2):59–60. doi: 10.1016/j.tree.2007.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

93. Khairallah M, Hoisington D. Laboratory protocols: CIMMYT applied molecular genetics laboratory. (1994). Ed. CIMMYT.

Google Scholar

94. Leelamanit W, Boonyom R, Panyim S, Reddy CC, Singh MM, Hayashi T, et al. A simple technique based on PCR to distinguish honeybee species. Bull Natl Institute Livestock Grassland Sci (Japan) (2003) 4:1–5.

Google Scholar

95. Beye M, Neumann P, Schmitzová J, Klaudiny J, Albert S, Šimúth J, et al. A simple, non-radioactive DNA fingerprinting method for identifying patrilines in honeybee colonies. Apidologie (1998) 29(3):255–63. doi: 10.1051/apido:19980305

CrossRef Full Text | Google Scholar

96. Hasselmann M, Fondrk MK, Page RE Jr., Beye M. Fine scale mapping in the sex locus region of the honey bee (Apis mellifera). Insect Mol Biol (2001) 10(6):605–8. doi: 10.1046/j.0962-1075.2001.00300.x

PubMed Abstract | CrossRef Full Text | Google Scholar

97. Lepais O, Weir J. SimRAD: an R package for simulation-based prediction of the number of loci expected in RADseq and similar genotyping by sequencing approaches an R package for simulation-based prediction of the number of loci expected in RADseq and similar genotyping by sequencing approaches. Mol Ecol Resources (2014) 14(6)1314–21. doi: 10.1111/1755-0998.12273

CrossRef Full Text | Google Scholar

98. Aguirre NC, Filippi CV, Zaina G, Rivas JG, Acuña CV, Villalba PV, et al. Optimizing ddRADseq in non-model species: A case study in Eucalyptus dunnii Maiden. Agronomy (2019) 9(9):484. doi: 10.3390/agronomy9090484

CrossRef Full Text | Google Scholar

99. Andrews S. FastQC: a quality control tool for high throughput sequence data (2010). Available at: http://www.bioinformatics.babraham.ac.uk/projects/fastqchttp://www.bioinformatics.babraham.ac.uk/projects/fastqc.

Google Scholar

100. Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA. Stacks: an analysis tool set for population genomics. Mol Ecol (2013) 22(11):3124–40. doi: 10.1111/mec.12354

PubMed Abstract | CrossRef Full Text | Google Scholar

101. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics (2014) 30(15):2114–20. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

102. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat met. (2012) 9(4):357–9. doi: 10.1038/nmeth.1923

CrossRef Full Text | Google Scholar

103. Yates AD, Allen J, Amode RM, Azov AG, Barba M, Becerra A, et al. Ensembl Genomes 2022: an expanding genome resource for non-vertebrates. Nucleic Acids Res (2022) 50(D1):D996–D1003. doi: 10.1093/nar/gkab1007

PubMed Abstract | CrossRef Full Text | Google Scholar

104. Yin L, Zhang H, Tang Z, Xu J, Yin D, Zhang Z, et al. rMVP: A memory-efficient, visualization-enhanced, and parallel-accelerated tool for genome-wide association study. Genom Proteom. Bioinform (2021) 19(4):619–28. doi: 10.1016/j.gpb.2020.10.007

CrossRef Full Text | Google Scholar

105. Zheng X, Levine D, Shen J, Gogarten SM, Laurie C, Weir BS. A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics (2012) 28(24):3326–8. doi: 10.1093/bioinformatics/bts606

PubMed Abstract | CrossRef Full Text | Google Scholar

106. R Core Team. R: a language and environment for statistical computing, version 3.0. 2. Vienna, Austria: R Foundation for Statistical Computing (2019) 2013.

Google Scholar

107. Ge SX, Jung D, Yao R. ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics (2020) 36(8):2628–9. doi: 10.1093/bioinformatics/btz931

PubMed Abstract | CrossRef Full Text | Google Scholar

108. Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H, et al. g: Profiler—a web server for functional interpretation of gene lists, (2016 update). Nucleic Acids Res (2016) 44(W1):W83–9. doi: 10.1093/nar/gkw199

PubMed Abstract | CrossRef Full Text | Google Scholar

109. Bu D, Luo H, Huo P, Wang Z, Zhang S, He Z, et al. KOBAS-i: intelligent prioritization and exploratory visualization of biological functions for gene enrichment analysis. Nucleic Acids Res (2021) 49(W1):W317–25. doi: 10.1093/nar/gkab447

PubMed Abstract | CrossRef Full Text | Google Scholar

110. Zheng X. SNPRelate: parallel computing toolset for genome-wide association studies. R Package version (2012) 95(10.18129):B9.

Google Scholar

111. Wickham H, Chang W, Wickham MH. Package ‘ggplot2’. Create elegant data visualisations using the grammar of graphics Vol. 2. United States: Springer-Verlag New York (2016). Available at: https://ggplot2.tidyverse.orghttps://ggplot2.tidyverse.org.

Google Scholar

112. Tallarida RJ, Murray RB, Tallarida RJ, Murray RB. Chi-square test. In: Manual of pharmacologic calculations: with computer programs, vol. 2012. Springer Science & Business Media (1987). p. 140–2.

Google Scholar

113. Mattila HR, Otis GW. Dwindling pollen resources trigger the transition to broodless populations of long-lived honeybees each autumn. Ecol Entomol (2007) 32(5):496–505. doi: 10.1111/j.1365-2311.2007.00904.x

CrossRef Full Text | Google Scholar

114. vanEngelsdorp D, Meixner MD. A historical review of managed honey bee populations in Europe and the United States and the factors that may affect them. J Invertebr. Pathol (2010) 103:S80–95. doi: 10.1016/j.jip.2009.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

115. Russo RM, Liendo MC, Landi L, Pietronave H, Merke J, Fain H, et al. Grooming behavior in naturally Varroa-resistant Apis mellifera colonies from north-central Argentina. Front Ecol Evol (2020) 8:590281. doi: 10.3389/fevo.2020.590281

CrossRef Full Text | Google Scholar

116. Otis GW. Revised distribution of three recently recognized species of honey bees in Asia. Honeybee Sci (1991) 15(4):167–70.

Google Scholar

117. Schneider SS, Lewis LA. The vibration signal, modulatory communication and the organization of labor in honey bees, Apis mellifera. Apidologie (2004) 35(2):117–31. doi: 10.1051/apido:2004006

CrossRef Full Text | Google Scholar

118. Guzmán-Novoa E, Hunt GJ, Page RE Jr., Fondrk MK. Genetic correlations among honey bee (Hymenoptera: Apidae) behavioral characteristics and wing length. Ann Entomol Soc Am (2002) 95(3):402–6. doi: 10.1603/0013-8746(2002)095[0402:GCAHBH]2.0.CO;2

CrossRef Full Text | Google Scholar

119. Withrow JM, Tarpy DR. Cryptic “royal” subfamilies in honey bee (Apis mellifera) colonies. PloS One (2018) 13(7):e0199124. doi: 10.1371/journal.pone.0199124

PubMed Abstract | CrossRef Full Text | Google Scholar

120. Brutscher LM, Baer B, Niño EL. Putative drone copulation factors regulating honey bee (Apis mellifera) queen reproduction and health: A review. Insects (2019) 10(1):8. doi: 10.3390/insects10010008

PubMed Abstract | CrossRef Full Text | Google Scholar

121. Lenoir JC, Laloi D, Dechaume-Moncharmont FX, Solignac M, Pham MH. Intra-colonial variation of the sting extension response in the honey bee Apis mellifera. Insectes Soc (2006) 53:80–5. doi: 10.1007/s00040-005-0838-5

CrossRef Full Text | Google Scholar

122. Zhang ZY, Li Z, Huang Q, Zhang XW, Ke L, Yan WY, et al. Deltamethrin impairs honeybees (Apis mellifera) dancing communication. Arch Environ Contam. Toxicol (2020) 78(1):117–23. doi: 10.1007/s00244-019-00680-3

PubMed Abstract | CrossRef Full Text | Google Scholar

123. Whitfield CW, Band MR, Bonaldo MF, Kumar CG, Liu L, Pardinas JR, et al. Annotated expressed sequence tags and cDNA microarrays for studies of brain and behavior in the honey bee. Genome Res (2002) 12(4):555–66. doi: 10.1101/gr.5302

PubMed Abstract | CrossRef Full Text | Google Scholar

124. Forêt S, Maleszka R. Function and evolution of a gene family encoding odorant binding-like proteins in a social insect, the honey bee (Apis mellifera). Genome Res (2006) 16(11):1404–13. doi: 10.1101/gr.5075706

PubMed Abstract | CrossRef Full Text | Google Scholar

125. Farooqui T, Vaessin H, Smith BH. Octopamine receptors in the honeybee (Apis mellifera) brain and their disruption by RNA-mediated interference. J Insect Physiol (2004) 50(8):701–13. doi: 10.1016/j.jinsphys.2004.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

126. Schlenstedt J, Balfanz S, Baumann A, Blenau W. Am5-HT7: molecular and pharmacological characterization of the first serotonin receptor of the honeybee (Apis mellifera). J Neurochem (2006) 98(6):1985–98. doi: 10.1111/j.1471-4159.2006.04012.x

PubMed Abstract | CrossRef Full Text | Google Scholar

127. Vleugels R, Verlinden H, Vanden Broeck J. Serotonin, serotonin receptors and their actions in insects. Neurotransmitter (2015) 2:e314. doi: 10.14800/nt.314

CrossRef Full Text | Google Scholar

128. Acevedo-Gonzalez JP, Galindo-Cardona A, Fuenzalida-Uribe NL, Ghezzi A, Giray T. Defensiveness measurement in honey bees (Apis mellifera) and brain expression of associated genes after noxious stimulus. bioRxiv (2022), 2022–04. doi: 10.1101/2022.04.15.488528

CrossRef Full Text | Google Scholar

129. Liu F, Shi T, Qi L, Su X, Wang D, Dong J, et al. lncRNA profile of Apis mellifera and its possible role in behavioural transition from nurses to foragers. BMC Genomics (2019) 20(1):1–11. doi: 10.1186/s12864-019-5664-7

PubMed Abstract | CrossRef Full Text | Google Scholar

130. Biswas S, Reinhard J, Oakeshott J, Russell R, Srinivasan MV, Claudianos C. Sensory regulation of neuroligins and neurexin I in the honeybee brain. PloS One (2010) 5(2):e9133. doi: 10.1371/journal.pone.0009133

PubMed Abstract | CrossRef Full Text | Google Scholar

131. Beggs KT, Hamilton IS, Kurshan PT, Mustard JA, Mercer AR. Characterization of a D2-like dopamine receptor (AmDOP3) in honey bee, Apis mellifera. Insect Biochem Mol Biol (2005) 35(8):873–82. doi: 10.1016/j.ibmb.2005.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

132. Wang Y, Jorda M, Jones PL, Maleszka R, Ling X, Robertson HM, et al. Functional CpG methylation system in a social insect. Science (2006) 314(5799):645–7. doi: 10.1126/science.1135213

PubMed Abstract | CrossRef Full Text | Google Scholar

133. Arakane Y, Dittmer NT, Tomoyasu Y, Kramer KJ, Muthukrishnan S, Beeman RW, et al. Identification, mRNA expression and functional analysis of several yellow family genes in Tribolium castaneum. Insect Biochem. Mol Biol (2010) 40(3):259–66. doi: 10.1016/j.ibmb.2010.01.012

CrossRef Full Text | Google Scholar

134. Paul RK, Takeuchi H, Kubo T. Expression of two ecdysteroid-regulated genes, Broad-Complex and E75, in the brain and ovary of the honeybee (Apis mellifera L.). Zool Sci (2006) 23(12):1085–92. doi: 10.2108/zsj.23.1085

CrossRef Full Text | Google Scholar

135. Dupuis JP, Gauthier M, Raymond-Delpech V. Expression patterns of nicotinic subunits α2, α7, α8, and β1 affect the kinetics and pharmacology of ACh-induced currents in adult bee olfactory neuropiles. J Neurophysiol (2011) 106(4):1604–13. doi: 10.1152/jn.00126.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

136. Ziai MR, Russek S, Wang HC, Beer B, Blume AJ. Mast cell degranulating peptide: a multi-functional neurotoxin. J Pharm- Pharmacol (1990) 42(7):457–61. doi: 10.1111/j.2042-7158.1990.tb06595.x

PubMed Abstract | CrossRef Full Text | Google Scholar

137. Labbé-Jullié C, Granier C, Albericio F, Defendini ML, Ceard B, Rochat H, et al. Binding and toxicity of apamin. Characterization of the active site. Eur J Biochem (1991) 196(3):639–45. doi: 10.1111/j.1432-1033.1991.tb15860.x

PubMed Abstract | CrossRef Full Text | Google Scholar

138. Welker S, Markert Y, Köditz J, Mansfeld J, Ulbrich-Hofmann R. Disulfide bonds of phospholipase A2 from bee venom yield discrete contributions to its conformational stability. Biochimie (2011) 93(2):195–201. doi: 10.1016/j.biochi.2010.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar