Fine-scale linkage mapping reveals a small set of candidate genes influencing honey bee grooming behavior in response to Varroa mites.

Populations of honey bees in North America have been experiencing high annual colony mortality for 15-20 years. Many apicultural researchers believe that introduced parasites called Varroa mites (V. destructor) are the most important factor in colony deaths. One important resistance mechanism that limits mite population growth in colonies is the ability of some lines of honey bees to groom mites from their bodies. To search for genes influencing this trait, we used an Illumina Bead Station genotyping array to determine the genotypes of several hundred worker bees at over a thousand single-nucleotide polymorphisms in a family that was apparently segregating for alleles influencing this behavior. Linkage analyses provided a genetic map with 1,313 markers anchored to genome sequence. Genotypes were analyzed for association with grooming behavior, measured as the time that individual bees took to initiate grooming after mites were placed on their thoraces. Quantitative-trait-locus interval mapping identified a single chromosomal region that was significant at the chromosome-wide level (p<0.05) on chromosome 5 with a LOD score of 2.72. The 95% confidence interval for quantitative trait locus location contained only 27 genes (honey bee official gene annotation set 2) including Atlastin, Ataxin and Neurexin-1 (AmNrx1), which have potential neurodevelopmental and behavioral effects. Atlastin and Ataxin homologs are associated with neurological diseases in humans. AmNrx1 codes for a presynaptic protein with many alternatively spliced isoforms. Neurexin-1 influences the growth, maintenance and maturation of synapses in the brain, as well as the type of receptors most prominent within synapses. Neurexin-1 has also been associated with autism spectrum disorder and schizophrenia in humans, and self-grooming behavior in mice.

Introduction

Ectoparasitic Varroa mites are considered by many to be the greatest threat to honey bee health worldwide. Increased annual mortality rates of North American colonies began about the time that tracheal mites (Acarapis woodii) and Varroa mites (V. destructor) first became established in the U.S. (about 1990) [1]. Varroa destructor switched hosts from Asian honey bees (Apis cerana) to the species commonly used for honey production and pollination worldwide (A. mellifera), about 60 years ago. Female Varroa mites lay their eggs within sealed brood cells in which bee larvae go through metamorphosis prior to emerging as adults. The mite progeny must mature and mate within the brood cell before the bee emerges for the mite to successfully reproduce. Virtually all honey bee colonies are infested with Varroa and unless steps are taken to reduce mite levels, colonies usually die within six months to two years, exhibiting dwindling populations and symptoms of viral and brood diseases [2].

Most colonies of bees in North America now show sufficient resistance to endoparasitic tracheal mites such that no treatment is needed to control them [3] but Varroa mite infestation has continued to be associated with high winter mortality in North America [1], [4]. The presence of V. destructor and the virus most commonly associated with mite infestation remain strong predictors of winter mortality in Europe [5]. Approximately 30% of bee colonies have been lost annually in the US in recent years [6]–[9]. Nationwide surveys of US bee losses were not published prior to 2006, but regional surveys taken after mites became established sometimes showed winter colony losses over 50%. Apiary losses were much reduced if beekeepers controlled Varroa mites [10]–[11]. As of 2009, high annual colony mortality has not been reported in Mexico [12]. However, Varroa is the health problem that has the greatest negative effect on honey production in Mexico and colony honey yields are significantly affected by mites [13]–[14]. Recently, surveys in North America and Europe also show a consensus that V. destructor infestation is still the most significant cause of annual colony mortality among all measured risk factors [15]–[20].

European races of A. mellifera are commonly used in beekeeping in most of the world, but an African race (A. mellifera scutellata) was introduced to Brazil in 1956. This race hybridized with European races, to form the Africanized honey bee, which has spread northward and has colonized all of Mexico and part of the US. The Africanized honey bee appears to have retained most of the original traits of the African race [21], including a somewhat higher degree of tolerance to Varroa than European races. But Africanized colonies still succumb to Varroa infestation [22]–[24]. A study conducted in Mexico identified grooming behavior as the most important factor that reduced mite population growth in a genetically diverse set of colonies. Colonies that had the lowest mite population growth during an eight-month period exhibited higher grooming behavior, had higher proportions of chewed mites falling from bees in the colonies, and reduced infestation levels of adult bees [25].

Mite-grooming behavior also varies between stocks of European honey bees in North America and can reduce infestation levels on adult bees [26]–[28]. The proportion of chewed mites falling from bees was found to correlate with removal rates in laboratory grooming assays [28], which corroborates earlier observations from hives in Europe [29]. Grooming behavior is an important resistance trait in Asian honey bees, the original host of Varroa, and is more strongly expressed in this species [30]–[31]. Recently, comparisons were made between four pairs of relatively mite-tolerant or susceptible honey bee stocks. The bees tested came from diverse geographic sources (Africanized bees tested in Mexico, Russian and European races tested in Canada). All mite-tolerant stocks showed comparitively higher proportions of chewed mites falling from colonies and increased intensity of individual grooming actions in laboratory assays, which underscores the importance of this trait [32].

Another honey bee mite-resistance mechanism was identified that involves bees uncapping infested brood cells, resulting in removal of pupae and mites from the cell or in the disruption of mite reproduction [33]–[35]. The latter trait has been called Varroa-sensitive hygiene, or VSH. Research has shown that Africanized honey bees exhibit a higher average levels of both grooming behavior and VSH compared to European races but it is unclear which trait is the most important for the Africanized honey bee’s increased tolerance to mites [22]–[24], [36].

Our objective was to use a QTL mapping approach to identifying candidate genes for honey bee mite-grooming behavior. In this report we describe fine-scale mapping of one putative QTL influencing mite-grooming behavior and the identification of a small set of candidate genes within the QTL region. Our approach was to use a single backcross family of worker bees derived from crosses between two colonies selected for high or low levels of grooming behavior. The time that individual beestook to exhibit grooming behavior after to contact with a single mite was quantified and these behavioral phenotypes were analyzed for association with alleles for about 1,300 intragenic single-nucleotide polymorphisms (SNPs). The results presented here and in a companion paper [37] are the first reports to identify candidate genes influencing behaviors that have been shown to be the key traits for suppressing mite population growth in screens of diverse populations of bees. These studies also have produced the most detailed and accurate linkage maps of the honey bee genome based on SNPs analyzed on genotyping arrays.

Results and Discussion

One putative QTL that we will refer to as groom-1 was identified with a LOD score of 2.72, covering about 2.0 Mb of chromosome 5 (Figure 1). On average, individuals that were homozygous for the high-grooming allele reacted to the presence of a mite faster than heterozygotes (18 versus 30 sec). Permutation tests indicated that this QTL exceeded the chromosome-wide significance level (p<0.05) but was not significant in a genome-wide test. Our results are relatively fine-scale in resolution; there were only 27 candidate genes that had been annotated in the official gene set 2 located in the LOD-1.5 support interval for this QTL, although future bioinformatic and EST analyses will probably identify additional genes and non-coding RNA transcripts in this region. The reason the candidate gene set is so small can be ascribed partly to the fact that this is not a gene-rich region, but is also the result of high marker saturation and the high recombination rate in the bee. The average recombination rate in the QTL region, 42 Kb/cM, was similar to the overall average for the bee, which has the highest reported rate for any metazoan species [38]–[39]. Two predicted genes could not be assigned any function (GB14792 and GB14110, Table 1). Two other genes, GB10440 and GB10743, are predicted to code for short proteins consisting of 80 and 242 amino acids, respectively, with high homology to proteins annotated in the bumble bee as ataxin-10 like. The presence of simple repeats within an intron of this gene in humans is associated with spinocerebellar ataxia type 10, characterized by atrophy of the cerebellum [40]–[41]. Two additional genes show homology to atlastin. The first is a putative gene designated GB15435, and annotated as containing only 81 amino acids. The other, GB14853, is the honey bee ortholog of atlastin-1, a membrane-bound large GTPase of the dynamin superfamily. Mutations in the human atlastin-1 are the most common risk factor for early onset of hereditary spastic paraplegia (HSP). Atlastin-1 protein is embedded in the endoplasmic reticulum (ER) and is critical for development of branched ER. In humans, defects in the ER are thought to be largely responsible for the etiology of HSP, which primarily affects long motor axons in humans that may extend up to several meters [42]. Atlastin-1 is one of several genes within the confidence interval that are important for function of the ER, microtubule function and protein trafficking (Table 1).

Figure 1

Map of QTL location for groom-1 on chromosome 5.

The physical location in base pairs of SNP probes in the honey bee genome assembly (Amel 4.0) is indicated to the right of the bar. The large number associated with the last SNP marker refers to its location in a contig that was not assigned to a chromosome prior to this study. Numbers to the left of the bar are distances in centimorgans. The red line indicates the LOD score for the likelihood that a QTL influencing grooming behavior is linked. The dotted line indicates the chromosome-wide empirical significance threshold of 0.05 as determined by 1000 permutations of phenotype data. Probe sequences matching the chromosome positions are available in .

Table 1

Candidate genes for honey bee mite-grooming behavior and their putative functions.

Gene name	Drosophila homolog	Predictions and proteindomains from blastp searches	Putative function
GB12154		Wnt-7b-like	Embryogenesis, morphogenic signaling
GB10440	CG4975	Ataxin-10-like	Neuron survival, differentiation, neurogenesis
GB10743	CG4975	Ataxin-10-like
GB16526	CG1093	Proteasome subunit alpha type-5 Ntn hydrolase	Fatty & amino acids metabolism
GB12705	CG30301	zinc finger protein 512B	Transcriptional regulation
GB15435		atlastin-like	Endoplasmic reticulum fusion, regulation of synaptic growth at neuromuscular junction
GB11239		Wnt-7b-like	Wnt signalling pathway
GB17462		regulator of microtubule dynamics protein 1-like	Retinol dehydrase 13-like
GB14853	CG6668	Atlastin-1	See GB15435
GB14672		ELM2 domain – unknown function SANT – DNA-binding domains	Transcriptional-regulating factor 1
GB18506	CG4406	Peptidase C13 superfamily	Putative GPI transamidase mediates GPI anchoring in the endoplasmic reticulum
GB13244		KAZAL_FS	Follistatin-like
GB10034	CG8079	G patch & FHA domains	angiogenic factor-like
GB17256	CG8014 Rme8	dnaJ homolog subfamily C member 13	Drosophila ortholog involved in receptor mediated endocytosis
GB20007		zinc finger protein 62-like THAP superfamily	THAP domain is a putative DNA-binding domain
GB17376		NADB Rossman superfamily FAR-C superfamily	fatty-acyl-CoA reductase 1-like
GB10140	CG12173	Phosphopentain phosphatase domain	Enolase-phosphatase E1; Methionine metabolism
GB19109		zinc finger protein 484-like	C2H2 domain, protein-protein interaction
GB16719	CG14181	Vesicle transport- USE1	Vesicle docking in the ER
GB17810	CG6176 Grip75	gamma tubulin complex component 4	Orientation of microtubules
GB18337	CG2774	sorting nexin 2-like	Intracellular protein trafficking
GB16547	CG18212	kinectin-like	Protein trafficking in the ER
GB11559	CG1972	Lactamase B superfamily RMMBL superfamily	RNA processing in translation and ribosome biogenesis (exonuclease)
GB18754	CG7050	Neurexin 1 EGF_CA and LNS superfamily domains	AmNrx1-Presynaptic protein involved in initiation, maintenance and function of synapses.
GB19804		Secapin preproprotein	Conserved venom peptide

The groom-1 QTL region also contains the sequence for the honey bee ortholog for neurexin-1, AmNrx1 (GB18754). As in other animals, AmNrx1 is highly expressed in the central nervous system and at a much lower level in other tissues. It is a large gene (∼400 Kb, 28 exons) that shows extensive alternative splicing. Twelve splice variants that differ in functional domains have been identified [43]. The two longest honey bee splice variants resemble human alpha-isoforms and one of the shortest forms resembles the beta isoform [44]. These two functional variants are generated through the use of alternative 3′ exons in the bee. The human gene generates these similar classes of isoforms by use of alternative promoters, suggesting convergent evolution. The intracellular domains of these two isoform classes are similar and in humans have been shown to interact with a number of proteins such as synaptogamin via PDZ domains, which can influence neurotransmitter release. The extracellular portion of the protein interacts with neuroligins embedded in the post-synaptic membrane and other synaptic proteins [45]–[49]. Nine isoforms in the bee lack the transmembrane domain. These apparently soluble isoforms may regulate neurexin functions and have been found in humans, although their functions apparently have not been determined [50]. Expression of AmNrx1 and neuroligin genes are concentrated in the mushroom body of the brain, which is the center of higher-order processing and learning in the bee [51]–[52]. AmNrx1 expression is also influenced by sensory experience, suggesting that it may play a role in the development of increased synaptic connections and behavioral plasticity [53].

The association of extracellular portions of neurexin and neuroligin proteins in nerve synapses serve important functions. These interactions form cellular adhesion complexes that regulate synapse formation, maintenance and maturation [45]–[48]. In humans and bees, the splice variants of neurexin-1 differ in the number of LNS motifs that interact with neuroligins in the synaptic gap. The longest alpha isoforms have three repeats, each consisting of an EGF (epidermal growth factor-like) sequence flanked by two LNS (laminin-neurexin-sex hormone binding globulin) motifs. The various isoforms of human neurexin-1 show specific differences in binding affinities for proteins in the synaptic cleft that regulate the function of the synapse [48]. Mutations of neurexin-1-alpha have been linked to autism spectrum disorders and schizophrenia in humans [54] and to synapse formation and associative learning in both Drosophila and the sea slug, Aplysia [55]–[56]. Knockout (KO) mice for neurexin-1-alpha were developed as a potential model system for autism spectrum disorder. Surprisingly, homozygous KO mice show greatly increased self-grooming activity compared to littermate controls, along with decreased nest building behavior and better performance in motor-learning tasks [57]. The grooming-behavior and motor-learning results suggest an increased sensitivity to tactile stimuli in the KO mice.

The possibility that neurexin-1 might influence grooming behavior in both honey bees and mammals is surprising, but of course discussion of candidate genes is speculative. This study identified just one putative QTL for honey bee mite-grooming behavior that was significant at the chromosome-wide level (LOD threshold = 2.28 for p<0.05) but not at the genome-wide level (LOD threshold = 3.49). It is necessary to confirm the effects of this QTL in a different population of bees. QTLs influencing behavioral traits usually have low LOD scores because these are traits influenced by large environmental effects and complex genetic architectures. In previous honey bee studies, only two QTL influencing the age of onset of foraging and one QTL influencing colony-level stinging response have met the threshold for genome-wide significance [58]–[59]. However, a number of honey bee behavioral QTL with low LOD scores have been confirmed in independent tests [60]. A recent study identified several QTL for a physiological mechanism of resistance in honey bees to Varroa mites but unfortunatelly these QTL only reached significance at the level of chromosomal region and the study could not narrow the list of candidate genes to a manageable number [61].

More effort is needed to improve the effectiveness of grooming assays. Our assay measured the time required for a bee to respond to the stimulus caused by the presence of a mite, which correlated with the proportion of mutilated mites falling from the test colonies (M. Arechavaleta-Velasco and K. Alcala-Escamilla, unpublished data). However, there was variation in how much the mites moved when on a bee. Perhaps an assay that measured the intensity of grooming actions in response to a biotic or abiotic stimulus might be better [32]. Differences in behavioral assays and/or seasonality in expression of the behavior are probably responsible for the discrepancy between heritability estimates and response to selection for mite-grooming traits [62]–[66]. It may also be useful to study correlated traits that may be influenced by the same causal gene(s). For example, one could test for effects on associative learning through training regimes that use the well-developed proboscis extension reflex assays. If the effect of the QTL on grooming behavior is confirmed, the identification of candidate genes that influence neural function in this region suggests promising directions for future research, such as correlating alternative splice variants or expression levels of AmNrx1 with honey bee mite-grooming behavior. RNAi knockouts of specific isoforms may also be possible if they can be targeted to the brain. If AmNrx1 does influence mite-grooming behavior, the honey bee would be a valuable comparative model to help understand behavioral and neurological influences of a protein that is highly conserved and involved in human neurological disorders. This research could also have an impact on breeding for resistance to mites.

Identification of a specific gene variant that has a major impact on behavioral resistance to Varroa mites would provide a means to more rapidly select for resistant strains of bees through a simple genotyping assay. It is important to actually identify the causal gene, or even the causal sequence variant, because the bee’s high recombination rate would likely result in many historical crossover events between a causal sequence and an intergenic DNA marker, reducing the linkage disequilibrium needed to predict the presence of favorable alleles. Marker-assisted breeding programs could be used to select colonies for high grooming behavior so that bees would be able to reduce their mite infestations to levels that would not affect honey production or compromise colony survival.

Materials and Methods

Ethics Statement

No permits were required to conduct the field research or genotyping analyses. The crosses and field research were conducted in Mexico under the supervision of researchers of the Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias, INIFAP Mexican agricultural research service). Genotyping was performed in the Purdue core genomics facility in accordance with university and federal biosafety regulations.

Behavioral Assay

Honey bee grooming behavior was measured in this study using the laboratory assay described by Espinoza et al. [66], with few modifications. The behavioral assay was conducted in an “observation frame” that was built using standard flat plastic hive frame foundation (45.0×16.0×1.5 cm) that has a pattern of hexagonal ridges and is coated with beeswax. This foundation served as the stage for observing bees. The foundation was inserted into a wooden hive frame of the appropriate size. The wood frame was covered on one of the ends with a piece of wire mesh of 16.0×10.0 cm for ventilation, the rest of the frame was covered with a transparent clear plastic sheet 1 mm thick. The wood frame was used to create a sealed area of 43.0×14.0×3.0 cm where bees could be confined an observed. In the center part of this area a circular area of 7.0 cm in diameter was built using wire mesh to create an inner arena inside the observation area, that excludes the passage of bees, but not the passage of air and to allow honeybee antennation and trophallaxis (touching of antennae and food sharing). A 2.5×2.5 cm window was cut in the plastic sheet above the circular area and a sliding piece of the plastic sheet was used as a cover to permit opening and closing of the window.

Adult Varroa mites were collected from a highly infested colony, by shaking adult bees from the hive frames into a 10 l plastic container, that had a 4.0 mm wire mesh screen 10 cm above its bottom. The container was closed with its plastic lid and CO2 was applied through a small hole in the cover and then the container was sealed for 3 min to anaesthetize the bees and to allow the mites to fall off from the workers. The mites were collected from the bottom of the container and placed in a petri dish and kept at room temperature in the laboratory. Approximately 300 worker bees from the colony to be tested were placed in the observation frame using a soft brush and four bees were placed in the inner circular arena. The observation frame was then sealed and taken to the laboratory. Sugar syrup (50% v/v) and water was fed to the bees through the wire mesh of the observation frame. An additional 40 honey bee workers were collected from the test colony, and each bee was introduced into a 10 ml Falcon™ tube and taken to the laboratory where they were kept at 8°C for a few minutes until their mobility was reduced to facilitate the manipulation of the bee.

The bee to be tested was transferred to the circular observation arena through the window in the plastic sheet with the help of a plastic tweezers. When the bee appeared to have fully recovered from the effect of the low temperature, a Varroa mite taken from the Petri dish was placed on her thorax through the window using a fine brush. The bee was observed by two people and the time that she took to respond by performing grooming behavior after the mite was placed over her body was recorded. Grooming behavior was defined as swiping motions in the direction of the mite with the front two pairs of legs. Colony phenotype was estimated by the mean reaction time of 40 workers.

Mapping Population

Two honey bee colonies were used as parental sources, one classified as a high groomer and one classified as a low groomer based on bi-directional selection in a population (n = 60) of honey bee colonies using the behavioral assay on 40 workers per hive. A queen was reared from each of the parental colonies and was artificially inseminated with the semen of three of her brothers. This inbreeding step was performed to ensure more genetically uniform F1 queens. A queen was reared from the high-grooming inbred colony and was artificially inseminated with the semen of a single haploid drone from the low-grooming inbred colony. From this queen, twenty F1 queens were reared and divided into two groups. Ten F1 queens were single-drone artificially inseminated with drones from the high grooming inbred colony and ten queens were single-drone inseminated with drones from the low-grooming inbred colony, to produce two types of colonies composed of backcross workers.

Each queen was introduced into a small colony made with three frames of brood, two frames of honey and approximately 1.5 kg of bees. The colonies were kept in single deep Jumbo type hives in the same apiary. All colonies were managed in the same way for a period of 60 days prior to the beginning of the experiments to allow time for workers in the colony to be replaced by daughters of the inseminated queens.

The grooming behavior of the 20 backcross colonies was tested using the behavioral assay, and the colony with the greatest range in grooming behavior phenotypes was selected for QTL mapping (a backcross to the high line). Worker bees (n = 400) from this colony were tested with the assay and the time that the each bee took to react by performing grooming behavior after a mite was placed on her thorax was recorded. At the end of each test the bee was collected in a 1.5 ml tube with ethanol and placed at −20°C. From the 400 bees that were tested, the 98 from the top of the phenotypic distribution and the 98 from the bottom of the distribution were selected for genotyping. Genomic DNA was extracted from each bee, DNA extraction involved grinding the bees in lysis solution (1% CTAB, 50 mM Tris pH 8.0, 10 mM EDTA, 1.1 M NaCl), followed by phenol-chloroform extraction and ethanol precipitation of the DNA [67]. The DNA of each individual bee was quantified and diluted in double distilled water to a final concentration of 50 ng/µl and kept at −80°C.

Genotyping and Linkage Analyses

Prior to genotyping the mapping population, the F1 queen that was the mother of the backcross colony was subjected to genomic sequencing on an ABI SOLiD platform. Heterozygous SNPs were identified using DiBayes software and these were compared to heterozygous SNPs in an unrelated F1 queen that was the mother of another mapping population (for a different QTL study) so that it was possible to design a set of probes that would be informative in both populations. Genotyping was performed using 250 ng of DNA of the 196 selected worker bees using the Illumina GoldenGate Assay, which usually achieves accuracies of ∼99.7 to 99.9% [68]. Briefly, DNA was fragmented and activated for binding to paramagnetic particles, then hybridized with allele-specific and locus-specific oligonucleotides. The last 3′nucleotide of the allele-specific nucleotide is at the SNP. Extension past the SNP and ligation to the locus-specific oligonucleotide follow, giving rise to full-length joined products that serve as templates for PCR with universal primers and dye-labeled allele-specific primers. This method allows for high levels of multiplexing. We used oligos to analyze 1,536 heterozygous SNPs in each individual. The dye-labeled PCR products were hybridized to the genotyping array matrix and the fluorescence signals were read by the BeadArray Reader and analyzed by Genome Studio software for semi-automated genotype clustering and calling.

SNP markers were assembled into linkage groups using JoinMap 4.0 software [69]–[70]. The marker orders were obtained by maximum likelihood analysis. Linkage distances between markers were estimated using the Kosambi mapping function. The time in seconds that it took for a bee to respond with grooming motions to a Varroa mite was log transformed to approximate a normal distribution. Interval mapping was performed with MapQTL 5.0 software [71]. Chromosome-wide and experiment-wide permutation tests in which phenotypes were randomly assigned to individuals were performed to calculate significance thresholds to identify significant and suggestive QTL. The 1.5 LOD support intervals, which approximately correspond to the 95% confidence intervals for the QTL position were determined from the interval mapping LOD values and the linkage map [72] and candidate genes in this genomic region were identified.

Analysis of Candidate Genes

The genes within the 1.5-LOD support interval were identified using the genome browser in BeeBase (http://hymenopteragenome.org/beebase/). The latest genome assembly scaffolds (Amel 4.5) within the 1.5-LOD interval were searched for annotated genes. Many of the genes are based on automated annotation software but considerable RNA sequence data has aided in these annotations. Blastp searches of the non-redundant protein databases provided a list of homologs and functional domains within the genes. Genes within the region were also assessed for putative function by reviewing the scientific literature concerning their homologs/orthologs.

Sequence of probes in .

(DOCX)

Click here for additional data file.