Skip to main content

Differential gene expression of the honey bee Apis mellifera associated with Varroa destructor infection



The parasitic mite, Varroa destructor, is the most serious pest of the western honey bee, Apis mellifera, and has caused the death of millions of colonies worldwide. This mite reproduces in brood cells and parasitizes immature and adult bees. We investigated whether Varroa infestation induces changes in Apis mellifera gene expression, and whether there are genotypic differences that affect gene expression relevant to the bee's tolerance, as first steps toward unravelling mechanisms of host response and differences in susceptibility to Varroa parasitism.


We explored the transcriptional response to mite parasitism in two genetic stocks of A. mellifera which differ in susceptibility to Varroa, comparing parasitized and non-parasitized full-sister pupae from both stocks. Bee expression profiles were analyzed using microarrays derived from honey bee ESTs whose annotation has recently been enhanced by results from the honey bee genome sequence. We measured differences in gene expression in two colonies of Varroa-susceptible and two colonies of Varroa-tolerant bees. We identified a set of 148 genes with significantly different patterns of expression: 32 varied with the presence of Varroa, 116 varied with bee genotype, and 2 with both. Varroa parasitism caused changes in the expression of genes related to embryonic development, cell metabolism and immunity. Bees tolerant to Varroa were mainly characterized by differences in the expression of genes regulating neuronal development, neuronal sensitivity and olfaction. Differences in olfaction and sensitivity to stimuli are two parameters that could, at least in part, account for bee tolerance to Varroa; differences in olfaction may be related to increased grooming and hygienic behavior, important behaviors known to be involved in Varroa tolerance.


These results suggest that differences in behavior, rather than in the immune system, underlie Varroa tolerance in honey bees, and give an indication of the specific physiological changes found in parasitized bees. They provide a first step toward better understanding molecular pathways involved in this important host-parasite relationship.


The honey bee (Apis mellifera, Insecta: Hymenoptera) has become an important model for genetic study, especially as its genome has been sequenced [1]. It is also an important economic insect as it is the world's principal crop pollinator and honey producer [2]. These activities have been threatened by the spread of Varroa destructor (Acari: Parasitiformes) a parasite of honey bees that causes devastating harm in many countries [3]. Varroa mites are ectoparasites of honey bees, parasitizing immature and adult bees and reproducing in cells in the honeycomb that contain brood [4]. Varroa mites or virus associated to mites impair the honey bee immune system [5] and in some cases boost the amplification of bee viruses [6]. The mechanisms underlying the mite's suppression of bee immunity and its impacts on pathogen virulence have not been elucidated.

There are genetic differences in the ability of honey bees to tolerate Varroa parasitism (tolerance is defined here as the capacity for the bees to survive when the parasite develops, preventing the demise of the hive). Some colonies of honey bees tolerate Varroa infection and survive despite the presence of the parasite in the hive [710]. These differences have been attributed to a variety of factors, including grooming and hygienic behavior [4], and differences in the timing of larval and pupal development that impact mite reproduction [11].

Colonies of the Asian honey bee Apis cerana (the original host of V. destructor[12]) suffer less damage from this parasite than A. mellifera in spite of the presence of the mite in the hives, and these two factors have been implicated [13]. The mechanisms underlying genetic differences for honey bee tolerance to Varroa mites are unknown. Insights into these mechanisms may lead to new molecular tools for both Varroa diagnosis and selective breeding of mite-tolerant honey bees for the bee industry. These issues are now amenable to study thanks to new genomic resources available for honey bees. Microarray analyses of differences in gene expression due to both mite parasitization and genotypic differences in bee tolerance are powerful approaches to explore the role of many genes in what are no doubt multifactorial resistance traits. Microarrays have proven to be useful in the study of host-pathogen interactions in other insects [1418].


Differences in bee gene expression in response to Varroa infestation and bee genotype

Comparisons of Varroa-infected and non-infected pupae from the 4 colonies studied pinpointed the expression of 32 genes whose expression varied significantly between the two types of pupae. Among them, 15 (47%) of these cDNAs were significantly up-regulated and 17 (53%) down-regulated in bees exposed to Varroa. The magnitude of the difference of expression is small in all cases. The single exception was EST (BB160020A20G03) which was over-expressed 20-fold in bees infected with Varroa. BLAST searches indicated that this EST matches the honey bee virus (AY292384), deformed wing virus [19].

There were 116 cDNAs whose expression varied with bee genotype, 47 (40.5%) of them up-regulated and 69 (59.5%) down-regulated. Two genes were regulated both as a function of Varroa infection and bee genotype. Dlic 2 is down-regulated in Varroa-parasitized bees and up-regulated in tolerant bees, whereas Strn-Mlck is down-regulated both in Varroa-parasitized and tolerant bees. The full list of genes with significant differences in expression is presented in Additional file 1.

Identification of A. mellifera-regulated transcripts and assignment of biological functions

To add biological meaning to the relatively large amount of microarray-derived data, we searched the GeneOntology database for the putative biological processes and molecular functions for the genes that showed significant differences in expression. Among these 148 genes, 19 associated with the presence of Varroa and 68 associated with bee genotype had information in Gene Ontology.

To unravel functional differences related to Varroa parasitism and bee genotype, we analyzed sets of genes for overrepresentation in functional GO categories according to biological process and molecular function (Additional file 2). Regarding biological process, genes involved in the regulation of protein metabolism, embryonic development and reproduction show a significant enrichment among the categories tested for the genes that varied in expression with Varroa infestation. A set of genes involved in cytokinesis, nervous system development, behavior, signal transduction and transcription-DNA dependent emerged from the analysis of genes that varied with bee genotype. The cluster analysis showing genes sharing the same GO term is shown in Table 1. Most of the genes affected by Varroa infection are involved in protein metabolism and their function is related to transferase and catalytic activity. The set of genes differentially expressed between tolerant and sensitive bees are mainly involved in transcription and neuron development. These genes have a nucleic acid binding function and are associated with pyrophosphatase GTPase, transferase and catalytic activity.

Table 1 Functional clustering based on Gene Ontology

RT-qPCR analysis

To further evaluate our ability to detect significant variation in gene expression at low fold-difference levels, a set of 4 functionally annotated bee transcripts was subjected to RT-qPCR analysis (Table 2). Differently expressed levels of cDNAs were quantified and are shown in Fig. 1. In accord with microarray analysis, cDNAs from Varroa-parasitized bees were up regulated for genes baz and CG9520 (unassigned) (Fig. 1a and 1b); or from tolerant bees for genes Alh and Hr78 (Fig. 1c and 1d). In all cases, expression levels quantified by qPCR are represented by the mean of 3 measurements (maximum and minimum values are indicated in Fig. 1).

Table 2 Characteristics of the RT-qPCR analysis of differentially expressed genes of Apis mellifera
Figure 1

cDNA levels of four transcripts shown by microarray analysis to be significantly regulated in Apis mellifera. In a and b, samples of bees infested with Varroa (hatched bars) are compared with bees free of Varroa (white bars). In c and d: Varroa-sensitive bees (hatched bars) are compared with Varroa-tolerant bees (white bars). All the values shown are mean ± SE.


Previous gene expression studies concerning honey bee immunity have mainly investigated responses to microbial pathogens [20]. Yet physiological responses to macroparasites such as Varroa are likely to be very different from microorganisms, as has been shown recently in a study on Drosophila revealing the ability of this species to activate a systemic immune response adapted to the invader [16]. The results described here pinpoint several genes regulated by Varroa parasitism that can be linked to honey bee responses to the presence of Varroa or to differences in bee tolerance.

Honey bee biological response to Varroa parasitism

The pathways suggested by this study to be associated with honey bee response to Varroa parasitism are schematized hypothetically in Fig. 2. Responses to Varroa can be grouped in two main categories, one related to deformed bee adults occasioned by the presence of the Deformed Wing Virus (DWV) often associated with Varroa parasitism, and a second category of genes that are related to cognitive impairment.

Figure 2

Hypothetical pathways and models of honey bee responses to Varroa -parasitism (A) and the bee tolerant genotype (B). Arrows and dashes indicate positive and negative regulation, respectively. Gene names in bold are up-regulated. In A) one of the consequences of the Varroa parasitism is a decline in immune capacity which induces the proliferation in bees of the Deformed Wing Virus (DWV). The boost of DWV multiplication might cause cellular and molecular damage, inducing the production of protein repair (Pcmt) and the labelling of proteins for degradation (Nedd8). In addition, regulated genes that might be affected by the presence of the DWV are indicated. Mites might decrease the production of dopamine (ple) and inhibit genes known for indirectly preventing neural degeneration in aged adults, which could explain the cognitive impairment often observed in adults parasitized by Varroa. In B) different genes can be associated to behavioral tolerance to Varroa. Tolerant and non-tolerant bees differ significantly by the expression of genes involved in the nervous system development. The olfactory pathway and neurons excitability seem also to play an important role in Varroa-tolerance. See text for a full discussion on the genes involved in the pathways presented.

Varroa parasitism and deformed wing virus

One of the consequences of Varroa parasitism is a decline in immune capacity which appears to induce the proliferation of viruses such as deformed wings virus in bees [6]. The down-regulation of the autophagic-specific gene 18 (Atg18) detected here is noticeable. Autophagy is an important mechanism for innate immunity against bacteria and viruses [21, 22]. The candidate innate immunity gene poly U binding factor 68 Kd (pUf68) [23] is also down-regulated in Varroa parasitized bees. By decreasing autophagy and immunity processes in bees, Varroa might favor the proliferation of DWV. Interestingly, the Varroa-parasitized bees displayed high levels of DWV viral RNA (about 20-fold; Additional file 1). The boost of DWV multiplication might cause cellular and molecular damage, and thus the observed production of genes for protein repair (Pcmt) [24], and the labeling of proteins for degradation (Nedd8) [25]. In contrast to these findings, we did not see a decrease on transcript abundance of immune pathway members [20] (Evans et al., 2006) found on this array. In fact, transcripts for the gene Rab7, a plausible regulator of immunity, were up-regulated in Varroa-parasitized bees.

The most notable symptoms of Varroa-parasitized bees are disfigured, small adults with deformed legs and wings [4, 6]. Our results show that at the transcriptome level the presence of Varroa down-regulates two genes involved in developmental processes, slg and dlg1. It has been shown that the Drosophila sugarless (sgl) gene regulates wingless signaling [26], which has a critical role in developmental processes. The gene dlg1 has been implicated in the control of proliferation of Drosophila imaginal discs programmed to produce adult structures at metamorphosis [27]. Although a down-regulation of these genes could have been induced by the presence of Varroa (via the DWV virus), a specific link with the development of deformed adults of virus-infected bees would need to be further investigated.

Cognitive impairment in bees parasitized by Varroa

Varroa infestation does not always cause wing deformity. Adults may sometimes appear to be normal morphologically, but there are mite effects on adult bee behavior. In particular, mite-parasitized foragers display a decrease in learning capability [28], prolonged absences from the nest and a lower rate of return to the colony [29]. A decrease in neuronal capacities involved in learning and navigation is a possible cause. Although the physiological mechanisms underlying reduced performance by bees in the presence of Varroa remain unknown, our results show that, the gene pale encoding tyramine hydroxylase is down-regulated in pupae parasitized by Varroa. Interestingly, the gene pale is needed for dopamine synthesis, which stimulates the nervous system and has many functions in the brain, including important roles in neural development, behavior and cognition, motor activity, motivation and learning. It is also interesting to note that the Dlic2 and Atg18 genes, both down-regulated in Varroa-parasitized bees, are enhancers of the blue cheese gene (bchs)[30], which is up-regulated in tolerant bees (see below). The bchs gene has been reported as preventing progressive neural degeneration in aged flies [31]. It is plausible, if these changes are chronic (still expressed in adults), that infested bees have a higher rate of neuronal apoptosis when aging, which might explain why foragers (the oldest bees in the nest) have difficulties in learning and orientating in flight.

Behavioral resistance of bees to Varroa

Grooming behavior seems to be one of the main characters involved in the tolerance of bees to Varroa [4]. Similarly, selection for tolerance to Varroa seems to be possible by selecting bees with a high level of hygienic behavior [32]. Hygienic behavior is understood as the ability of bees to uncap and remove infected brood. Although the molecular mechanisms underlying this trait remain unknown, some insights are provided here based on the large-scale comparative analysis of transcriptome differences between Varroa-tolerant and sensitive bees.

A disproportionately high fraction of the genes differentially expressed between tolerant and susceptible bees are involved in the development of the nervous system (Fig. 2B). A large part of these genes are down-regulated in tolerant bees compared to sensitive bees, including the genes futsch, scratch (scrt), otk, Myosin heavy-chain-like (Mhcl), groucho (gro), kekkon-1 (kek1) [26, 3338]. Fringe (fng), also down regulated in Varroa-tolerant bees, plays an important role in the positive regulation of Notch signalling pathway involved in the regulation of genes that control multiple cell differentiation processes during embryonic and adult life, such as neuronal function and development [39]. Finally, a gene involved in locomotory behavior, single-minded (sim), is down-regulated in tolerant bees. Drosophila mutant sim flies are only able to walk in circles and this phenotype is due to defects in the central brain complex [40].

Several genes involved in neuron excitability are up-regulated in Varroa-tolerant bees compared to sensitive bees. Purity of essence (poe) (also named pushover) is involved in neuronal excitability and may play a role in responsiveness to environmental stimuli and behavior. For example, fly mutants display behavioral defects such as sluggishness, uncoordination and defective flight [41]. This gene has been shown to be down-regulated in honey bees by the queen mandibular pheromone, which slows the behavioral maturation of workers, e.g. transition from inside hive work (nurse) to foraging activity(forager) [42]. The GluClα gene, a glutamate-gated Cl- channel specific to arthropods, which is also up-expressed in tolerant bees, is known to modulate neuronal membrane excitability [43]. The paralytic gene (para), known to be important in the conducting of nerve action potentials in flies [44], is also up-regulated. Seen together, these genes suggest a mechanism by which Varroa-tolerant bees are more sensitive to external stimuli than Varroa sensitive bees. Also, the Dynein heavy chain 64C gene (Dhc64c), up-regulated in tolerant bees, is required for proliferation of mushroom-body neuroblasts [45]. Mushroom bodies have an important role in insect cognition and are known to be involved in learning and memory, particularly for smell.

Interestingly, several genes involved in olfaction (smell impaired 21F smi21F, suppressor of the white-apricot gene su(wa), poe, para, Rogdi) are up-regulated in Varroa-tolerant bees compared to sensitive bees [4650]. However, the Down syndrome cell adhesion molecule gene (Dscam) and four wheel drive (fwd) also have a role in olfaction, and they are down-regulated in Varroa-tolerant bees [34]. The differential expression of this group of genes depending on bee genotype is of major importance considering that in the hive, olfaction and neuronal sensitivity together may play a major role in the detection of Varroa infested cells. Observations made on two earlier generations of the bee stocks here studied have been shown that Varroa-tolerant bees have better ability to detect the mite [51]. Previous studies showed that hygienic bees have a higher olfactory sensitivity and responsiveness than non-hygienic bees [5254]. They are notably able to discriminate between odors of healthy and diseased brood at a lower stimulus level, suggesting that olfaction and responsiveness play a key role in hygienic behavior. In addition, observations made on the behavior of the brood, have shown open cells containing destroyed parasitized bee pupae in tolerant experimental colonies, which is in agreement with a bee hygienic behavior against Varroa. If these gene candidate pathways for Varroa behavioral tolerance can be further confirmed, the results suggest that they already set up during the pupae stage.

There were differences in expression between tolerant and sensitive bees for genes involved in increased resistance to toxins like Ahcy13 (Adenosylhomocysteinase at 13) involved in detoxification [55], and para involved in resistance to insecticide pyrethroids and DDT (Dichloro-Diphenyl-Trichloroethane) [56]. Two genes linked to immunity, Dsam and otk, members of the Immunoglobulin gene superfamily in Drosophila [57] were down-regulated in Varroa-tolerant bees, while one galectin-family gene (an apparent ortholog to Drosophila CG32226) appears to be upregulated.

Low fold differences in gene expression

For most of the genes showing significant differences in expression in the present study, the magnitude of these differences was small. Recent reports have shown that microarrays can significantly underestimate gene expression changes and therefore if a severe cut-off is applied, this approach might miss important changes in gene expression. Recent reports have validated the capability of a microarray approach to detect small gene expression differences [58]; 87% of a set of parasite-specific genes displayed changes 2-fold or less in D. melanogaster challenged to a protozoan parasite [17]. Similarly, about 60% of genes of the lepidopteran, Spodoptera frugiperda, revealed around 0.5-fold changes in the transcript levels associated with virus infection [14]. As one possible explanation for low level gene regulation, it has been suggested in Drosophila that much of the response to parasite attack probably does not involve de novo gene expression but post-transcriptional events [18]. Another possibility is that many small differences in gene expression reflect subtle modulation by a large number of factors acting in cascade. The small changes reported here might also represent an underestimation of tissue-specific effects obscured by whole-body analysis. Future studies of Varroa effects on honey bees should analyze specific tissues, and our results suggest a focus on the brain would be fruitful. Previous studies have demonstrated extensive regulation of brain gene expression in conjunction with honey bee behaviour [42, 5961].


This work is the first step towards understanding the genomic responses of honey bees to Varroa parasitism. It demonstrates that honey bee pupae exhibit differences in gene expression associated either with the presence of Varroa, or with tolerance to this parasite. These results highlight the potential importance of behavioral mechanisms of response to Varroa and suggest that a study focused on the brain is of importance for the future. For an economically important species such as the honey bee, the identification of parasite-specific response factors might ultimately serve to identify molecules that act on bee parasites. In addition, differences between tolerant and sensitive bees could lead to developing tools to select improved strains of honey bees for beekeepers.


Honey bee colonies and sample collection

Both, sensitive and tolerant honey bee colonies here studied, belong to non related local strains of the same A. mellifera population which is bred in the Laboratory of Bee Biology and Protection, Institut National de la Recherche Agronomique, Avignon, France. Four colonies were used for this study; two (T222 and T757) had survived for 11 years with no chemical treatment against Varroa despite the widespread presence of this parasite in the locale, and displayed a very low rate of parasitism [10]. The other two colonies (S21 and S38) showed a very high level of parasitism. No acaricide treatment was applied to avoid pesticide bias. Parasitism intensity was determined by counting the number of mites that died naturally and accumulated on the floor of the hive in each colony during the year (method described in Ellis [62]). Estimations based on data for 5 years (4 measurements/year) showed that the infestation rate in colonies S21 and S38 was on average 10 times higher than in colonies T222 and T757.

The reproduction of Varroa mites can be affected by traits in developing bees [63, 64]. We therefore examined gene expression in pupae. Honey bee pupae were collected from the four bee colonies at the blue-eye stage at the start of cuticle pigmentation, making it possible to determine the bee developmental stage [65]. Capped brood cells were opened and two samples of 50 parasitized and non-parasitized pupae were collected. A pupa was considered to be parasitized if a reproductive Varroa was found with it in the honeycomb cell. The pupae were collected and snap frozen using N2 and stored at -80°C until RNA extraction.

We compared parasitized and non-parasitized pupae from both susceptible and tolerant colonies. To minimize the effect of genetic variation, pools of full-sister pupae (related by 75% due to haplodiploidy [66]) were collected in equal numbers from each of the study colonies.

Sibling estimation by microsatellite genotyping

To identify full sisters, one hundred bees per colony were genotyped (50 parasitized and 50 not parasitized) using Ap53 and A107 microsatellite loci; these loci have been shown to be effective for this purpose [67]. DNA was extracted from 2 legs/pupa, and DNA extraction and PCR amplification were as described [67], except that the PCR products were detected on a MegaBACE DNA Analysis System 1000 (Molecular Dynamics Inc., USA). The genotype data were used to assemble sets of 8 pupae belonging to the same full sibling group. For each of these, subsets of pupae of Varroa-parasitized and non parasitized bees were analyzed in microarray experiments. This procedure minimizes the effects of intra-colonial genetic variation on gene expression.

Microarray experimental design

We compared parasitized and non-parasitized pupae from both susceptible and tolerant colonies using a direct loop design. Direct comparisons were made between parasitized and unparasitized bees from sensitive colonies and from tolerant colonies (Fig. 3). All comparisons were performed in dye-swap on 2 biological samples. A total of 8 samples were thus contrasted using 16 arrays in all (2 biological replicates and 2 technical replicates). Each sample consisted of a fraction of the total RNA obtained from a pool of 8 bees. Full-sister bees were assembled as pools from specific patrilines detected in each colony.

Figure 3

Microarray experimental design. For all experiments, arrows indicate microarray hybridizations (arrow tail, Cy3-labeled sample; arrowhead, Cy5-labeled sample). Varroa parasitized (+) and non-parasitized (-) full-sister pupae, from two different genetic stocks: one susceptible (S) and one tolerant (T) to Varroa were compared. Dye swaps were made for all comparisons.

RNA isolation

The sets of bees were freeze-dried and ground in liquid nitrogen. 30 mg (about 3% of the powder obtained) was used for each RNA extraction using the RNeasy Plant Mini Kit (Qiagen, Courtaboeuf, France) following the manufacturer's protocol for animal tissues. The total RNA solution was treated in liquid by RNase-Free DNase I and purified on column RNA Cleanup (Qiagen, Courtaboeuf, France) according to the manufacturer's protocols. The purity and concentration of RNA were determined by OD measurements in a spectrophotometer. RNA extraction was validated by specific PCR amplifications following RT-PCR of both the glyceraldehyde 3 phosphate dehydrogenase 1 (XM_397363) and the epsilon-tubulin 1 (XM_394700) genes. The quality of the extracted RNA (integrity and size distribution of total RNA) was verified by 1.2% agarose gel electrophoresis in TAE buffer and ethidium bromide staining.

cDNA synthesis and array hybridizations

10 μg total RNA was used for cDNA synthesis using fluorescent Cy3-dCTP or Cy5-dCTP (Amersham Biosciences, Orsay, France) and the Pronto!TM plus System kit (Promega, Charbonnières, France) for labeling. Dye-labeled cDNAs were purified on cleanup columns (Promega, Charbonnières, France). 55 pmol of the labeled cDNA per dye was used per slide. Prehybridation, hybridization and washing steps were performed according to [61]. The slides were hybridized at 42°C and air dried by centrifugation for 2 min at 800 rpm at room temperature.

Microarray data acquisition and statistical analysis

The arrays were scanned on a GenePix 4000A scanner (Axon Instruments, Foster City, USA) and images were analyzed by GenePix Pro 3.0 (Axon Instruments, Foster City, USA). For each array, the raw data comprised the logarithm (base 2) of median feature pixel intensity at wavelength 635 nm (red) and 532 nm (green). No background was subtracted. We excluded control spots (as described in [61]) and spots for which duplicates did not pass quality controls standards, cDNAs were not found by spot-finding software, or those determined to be irregular by visual inspection. Intensity signals for cDNAs passing these filters were normalized for intensity- and position-dependent bias. Array-by-array normalization was performed to remove systematic biases. Then, we replaced the value of the spots that were considered as badly formed features with the value of the duplicate. We averaged the two values from each duplicated feature to obtain one value of the gene per array in each condition.

A total of 4,795 cDNAs passed initial filters; 3,045 of these were collapsed to genes and 1,750 ESTs were unassigned. We refer to the combined set of genes and unassigned ESTs as genes, although some redundancy might exist. ANOVA analysis was used to classify the differentially expressed genes according with the different factors considered (P < 0.05 was used to denote statistical significance).

Statistical analyses were conducted using R version 2.2.0 [61] and the R/maanova software package version 0.98.8 [61, 68]. The gene-by-gene mixed effect ANOVA model, Yijkvr = μ + Ai + Dj + (ST)k + Vv + Eijkvr, was applied. Observation Yijkvr is the expression value of the gene in the biological replicate r studied on array i when the condition is labelled with dye j. The condition is defined by a modality of the genotype, say k, and the modality of the Varroa, say v. The residual of the mode Eijkvr and the term Ai are treated as random effects and the others as fixed effects. Statistical tests were performed with R/Maanova using the hybrid variance model, Fs, based on a James-Stein estimator [69, 70]. To control the number of false-positives, the significance (nominal P-value) of genes was computed using a permutation test (pvalperm with 1000 permutations) [70]. Permutation tests have the advantage of not assuming a parametric underlying distribution of expression values. A gene was declared differentially expressed if its adjusted P-value is lower than 0.05. The adjusted P-value used here made it possible to control the Family Wise Error Rate (FWER)

Functional analysis of gene expression

Microarray ESTs were annotated as described in [60, 71]. Briefly, ESTs corresponding to microarray cDNAs were tested for near-perfect matches (98% identity) to coding (protein) sequence or to genomic sequence within or immediately downstream (500 bp) of predicted genes (using release 1 of the honey bee Official Gene Set [72]). Redundant cDNA values were averaged (by using untransformedvalues), and resulting values were assigned to official gene names (which are all prefixed "GB"). Remaining cDNAs not associated with predicted sequences retained their EST identifiers and are presented here by EST accession number (prefixed "BB"). Genes were tentatively assigned molecular function terms based on annotation of the single best BLASTX match in Drosophila melanogaster.

We used Gene Ontology (GO) analysis to identify the biological function of the differentially expressed genes, as described in [60]. We sought statistically overrepresented terms among the set of genes significantly regulated. This Enrichment of Gene Ontology (GO) category was tested with a Hypergeometric test followed by the Benjamini Hochberg correction for multiple testing using the analysis tool GOToolBox [60]. Briefly, the analysis calculates the frequency of terms of a gene list (the significantly regulated genes) and compares these results with total frequencies of genes analyzed. We also grouped functionally related genes on the basis of their GO terms. Distances based on GO terms are calculated for all possible pairs from the gene list which are then used for clustering (e.g. the probability that they are functionally related). The functional clustering of genes was carried out in GOToolBox using the WPGMA algorithm. A minimum cluster size of 3 genes was applied.

Transcript quantification by RT-qPCR

The results of the microarray study were confirmed by measuring the expression of 4 differentially regulated genes (Table 2) using Real-Time quantitative PCR (RT-qPCR). Using an aliquot of the RNA extractions (pools of 8 bees) used for the microarray work, cDNAs were created using the ImProm-II™ Reverse Transcription System kit (Promega). This cDNA (three replicates) was used as template for the Real-Time PCR performed in a Roche LightCycler 480 Instrument platform. The reaction mix consisted of 10 μl of LightCycler 480 SYBR Green Supermix (Roche Laboratories), 5 μl (either 600 nM or 300 nM) of forward and reverse primers, 3 μl dH2O and 2 μl cDNA template. Three reactions were performed and means calculated for each locus and treatment. To standardize the results, a housekeeping gene (Rp49) that did not vary in expression was used as control. The forward/reverse primer sequences are indicated in Table 2. The efficiency for each locus was determined by running a dilution series (1000x, 100x, 10x, 1x) in triplicate. The results were standardized using the [60] method. Efficiency of the amplicons obtained for each locus was adequately high and at least 97% (Table 2).


  1. 1.

    Consortium HGS: Insights into social insects from the genome of the honeybee Apis mellifera. Nature. 2006, 443 (7114): 931-949.

    Article  Google Scholar 

  2. 2.

    Burgett M, Rucker RR, Thurman WN: Economics and honey bee pollination markets. Am Bee J. 2004, 144 (4): 269-271.

    Google Scholar 

  3. 3.

    Griffiths DA, Bowman CE: World distribution of the mite Varroa jacobsoni, a parasite of honeybees. Bee World. 1981, 62: 154-163.

    Article  Google Scholar 

  4. 4.

    Sammataro D, Gerson U, Needham G: Parasitic mites of honey bees: life history, implications and impact. Annu Rev Entomol. 2000, 45: 519-548.

    PubMed  CAS  Article  Google Scholar 

  5. 5.

    Gregory PG, Evans JD, Rinderer TE, De Gruzman L: Conditional immune-gene suppression of honeybees parasitized by Varroa mites. J Insect Scien. 2005, 5: 1-5.

    Article  Google Scholar 

  6. 6.

    Yang XY, Cox-Foster DL: Impact of an ectoparasite on the immunity and pathology of an invertebrate: evidence for host immunosuppression and viral amplification. Proc Natl Acad Sci USA. 2005, 102 (11): 7470-7475.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  7. 7.

    Fries I, Camazine S, Sneyd J: Population-Dynamics of Varroa jacobsoni – a Model and a Review. Bee World. 1994, 75 (1): 5-28.

    Article  Google Scholar 

  8. 8.

    Seeley TD: Honey bees of the Arnot Forest: a population of feral colonies persisting with Varroa destructor in the northeastern United States. Apidologie. 2007, 38 (1): 19-29.

    Article  Google Scholar 

  9. 9.

    Wood M, Hardin B, Lee J: Varroa-tolerant bees keep hives buzzing. Agr Res Mag. 1999, 47: 10-12.

    Google Scholar 

  10. 10.

    Le Conte Y, de Vaublanc G, Crauser D, Jeanne F, Rousselle JC, Bécard JM: Honey bee colonies that have survived Varroa destructor. Apidologie. 2007, 38:

    Google Scholar 

  11. 11.

    Moritz FA: Heritability of the postcapping stage in Apis mellifera and its relation to varroatosis resistance. The Journal of Heredity. 1985, 76: 267-270.

    Google Scholar 

  12. 12.

    Anderson DL, Trueman JWH: Varroa jacobsoni (Acari: Varroidae) is more than one species. Exp Appl Acarol. 2000, 24 (3): 165-189.

    PubMed  CAS  Article  Google Scholar 

  13. 13.

    Peng YS, Y F, S X, L G: The resistance Mecanism of the Asian Honey Bee, Apis cerana Fabr., to an Ectoparasitic Mite, Varroa jacobsoni Oudemans. J Invert Path. 1987, 49: 54-60.

    Article  Google Scholar 

  14. 14.

    Barat-Houari M, Hilliou F, Jousset FX, Sofer L, Deleury E, Rocher J, Ravallec M, Galibert L, Delobel P, Feyereisen R, Fournier P, Volkoff AN: Gene expression profiling of Spodoptera frugiperda hemocytes and fat body using cDNA microarray reveals polydnavirus-associated variations in lepidopteran host genes transcript levels. BMC Genomics. 2006, 7: (Jun. 21): Art. No. 160.

    Google Scholar 

  15. 15.

    Sim C, Hong YS, Vanlandingham DL, Harker BW, Christophides GK, Kafatos FC, Higgs S, Collins FH: Modulation of Anopheles gambiae gene expression in response to o'nyong-nyong virus infection. Insect Mol Biol. 2005, 14 (5): 475-481.

    PubMed  CAS  Article  Google Scholar 

  16. 16.

    Vodovar N, Vinals M, Liehl P, Basset A, Degrouard J, Spellman P, Boccard F, Lemaitre B: Drosophila host defense after oral infection by an entomopathogenic Pseudomonas species. Proc Natl Acad Sci USA. 2005, 102 (32): 11414-11419.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  17. 17.

    Roxstrom-Lindquist K, Terenius O, Faye I: Parasite-specific immune response in adult Drosophila melanogaster: a genomic study. Embo Reports. 2004, 5 (2): 207-212.

    PubMed  PubMed Central  Article  Google Scholar 

  18. 18.

    Wertheim B, Kraaijeveld AR, Schuster E, Blanc E, Hopkins M, Pletcher SD, Strand MR, Partridge L, Godfray HCJ: Genome-wide gene expression in response to parasitoid attack in Drosophila. Genome Biology. 2005, 6 (11): Art. No. R94.

    Google Scholar 

  19. 19.

    Lanzi G, Miranda JR, Boniotti MB, Cameron CE, Lavazza A: Molecular and biological characterization of deformed wing virus of honeybees (Apis mellifera L.). J Virol. 2006, 80: 4998-5009.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  20. 20.

    Evans JD: Beepath: an ordered quantitative-PCR array for exploring honey bee immunity and disease. J Invertebr Pathol. 2006, 93 (2): 135-139.

    PubMed  CAS  Article  Google Scholar 

  21. 21.

    Deretic V: Autophagy in innate and adaptive immunity. Trends in Immunology. 2005, 26 (10): 523-528.

    PubMed  CAS  Article  Google Scholar 

  22. 22.

    Levine B: Eating oneself and uninvited guests: autophagy-related pathways in cellular defense. Cell. 2005, 120 (2): 159-162.

    PubMed  CAS  Google Scholar 

  23. 23.

    Foley E, O'Farrell PH: Functional dissection of an innate immune response by a genome-wide RNAi screen. PLoS Biology. 2004, 2 (8): 1091-1106.

    CAS  Article  Google Scholar 

  24. 24.

    O'Connor MB, Galus A, Hartenstine M, Magee M, Jackson FR, O'Connor CM: Structural organization and developmental expression of the protein isoaspartyl methyltransferase gene from Drosophila melanogaster. Insect Biochem Mol Biol. 1997, 27 (1): 49-54.

    PubMed  Article  Google Scholar 

  25. 25.

    Ou C, Lin Y, Chen Y, Chien C: Distinct protein degradation mechanisms mediated by Cul1 and Cul3 controlling Ci stability in Drosophila eye development. Genes Dev. 2002, 16: 2403-2414.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  26. 26.

    Hacker U, Lin X, Perrimon N: The Drosophila sugarless gene modulates Wingless signaling and encodes an enzyme involved in polysaccharide biosynthesis. Development. 1997, 124 (18): 3565-3573.

    PubMed  CAS  Google Scholar 

  27. 27.

    Woods D, Hough C, Peel D, Callaini G, Bryant P: Dlg protein is required for junction structure, cell polarity, and proliferation control in Drosophila epithelia. J Cell Biol. 1996, 134: 1469-1482.

    PubMed  CAS  Article  Google Scholar 

  28. 28.

    Kralj J, Brockmann A, Fuchs S, Tautz J: The parasite mite Varroa destructor affects non-associative learning in honey bee foragers, Apis mellifera L. J Comp Physiol A Neuroethol Sens Neural Behav Physiol. 2007, 193 (3): 363-370.

    PubMed  Article  Google Scholar 

  29. 29.

    Kralj J, Fuchs S: Parasitic Varroa destructor mites influence flight duration and homing ability of infested Apis mellifera foragers. Apidologie. 2006, 37 (5): 577-587.

    Article  Google Scholar 

  30. 30.

    Simonsen A, Cumming RC, Lindmo K, Galaviz V, Cheng S, Rusten TE, Finley KD: Genetic modifiers of the Drosophila blue cheese gene link defects in lysosomal transport with decreased life span and altered ubiquitinated-protein profiles. Genetics. 2007, 176 (2): 1283-1297.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  31. 31.

    Finley KD, Edeen PT, Cumming RC, Mardahl-Dumesnil MD, Taylor BJ, Rodriguez MH, Hwang CE, Benedetti M, McKeown M: Blue cheese mutations define a novel, conserved gene involved in progressive neural degeneration. Journal of Neuroscience. 2003, 23 (4): 1254-1264.

    PubMed  CAS  PubMed Central  Google Scholar 

  32. 32.

    Harbo JR, Harris JW: Suppressed mite reproduction explained by the behavior of adult bees. J Apic Res. 2005, 44 (1): 21-23.

    Google Scholar 

  33. 33.

    Cafferty P, Yu L, Rao Y: The receptor tyrosine kinase off-track is required for layer-specific neuronal connectivity in Drosophila. Development. 2004, 131 (21): 5287-5295.

    PubMed  CAS  Article  Google Scholar 

  34. 34.

    Hummel T, Vasconcelos M, Clemens J, Fishilevich Y, Vosshall L, Zipursky S, Anholt ea: Axonal Targeting of Olfactory Receptor Neurons in Drosophila Is Controlled by Dscam. Neuron. 2001, 37: 221-231.

    Article  Google Scholar 

  35. 35.

    Ivanov AI, Rovescalli AC, Pozzi P, Yoo S, Mozer B, Li HP, Yu SH, Higashida H, Guo V, Spencer M, Nirenberg M: Genes required for Drosophila nervous system development identified by RNA interference. Proc Natl Acad Sci USA. 2004, 101 (46): 16216-16221.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  36. 36.

    Musacchio M, Perrimon N: The Drosophila kekkon genes: novel members of both the leucine-rich repeat and immunoglobulin superfamilies expressed in the CNS. Dev Biol. 1996, 178 (1): 63-76.

    PubMed  CAS  Article  Google Scholar 

  37. 37.

    Paroush Z, Finley RL, Kidd T, Wainwright SM, Ingham PW, Brent R, Ish-Horowicz D: Groucho is required for Drosophila neurogenesis, segmentation, and sex determination and interacts directly with hairy-related bHLH proteins. Cell. 1994, 79: 805-815.

    PubMed  CAS  Article  Google Scholar 

  38. 38.

    Roark M, Sturtevant MA, Emery J, Vaessin H, Grell E, Bier E: Scratch, a pan-neural gene encoding a zinc finger protein related to snail, promotes neuronal development. Genes & Development. 1995, 9 (19): 2384-2398.

    CAS  Article  Google Scholar 

  39. 39.

    Bolos V, Grego-Bessa J, de la Pompa JL: Notch signaling in development and cancer. Endocr Rev. 2007, 28 (3): 339-363.

    PubMed  CAS  Article  Google Scholar 

  40. 40.

    Pielage J, Steffes G, Lau DC, Parente BA, Crews ST, Strauss R, Klambt C: Novel behavioral and developmental defects associated with Drosophila single-minded. Dev Biol. 2002, 249 (2): 283-299.

    PubMed  CAS  Article  Google Scholar 

  41. 41.

    Richards S, Hillman T, Stern M: Mutations in the Drosophila pushover gene confer increased neuronal excitability and spontaneous synaptic vesicle fusion. Genetics. 1996, 142 (4): 1215-1223.

    PubMed  CAS  PubMed Central  Google Scholar 

  42. 42.

    Grozinger CM, Sharabash NM, Whitfield CW, Robinson GE: Pheromone-mediated gene expression in the honey bee brain. Proc Natl Acad Sci USA. 2003, 100 (Suppl 2): 14519-14525.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  43. 43.

    Palladino M, Keegan L, O'Connell M, Reenan R: A-to-I Pre-mRNA Editing in Drosophila Is Primarily Involved in Adult Nervous System Function and Integrity. Cell. 2000, 102: 437-449.

    PubMed  CAS  Article  Google Scholar 

  44. 44.

    Ganetzky B: Genetic studies of membrane excitability in Drosophila: lethal interaction between two temperature-sensitive paralytic mutations. Genetics. 1984, 108 (4): 897-911.

    PubMed  CAS  PubMed Central  Google Scholar 

  45. 45.

    Liu Z, Steward R, Luo L: Drosophila Lis1 is required for neuroblast proliferation, dendritic elaboration and axonal transport. Nature Cell Biology. 2000, 2: 776-783.

    PubMed  CAS  Article  Google Scholar 

  46. 46.

    Anholt RR, Fanara JJ, Fedorowicz GM, Ganguly I, Kulkarni NH, Mackay TF, Rollmann SM: Functional genomics of odor-guided behavior in Drosophila melanogaster. Chem Senses. 2001, 26 (2): 215-221.

    PubMed  CAS  Article  Google Scholar 

  47. 47.

    Anholt RR, Lyman RF, Mackay TF: Effects of single P-element insertions on olfactory behavior in Drosophila melanogaster. Genetics. 1996, 143 (1): 293-301.

    PubMed  CAS  PubMed Central  Google Scholar 

  48. 48.

    Dubnau J, Chiang AS, Grady L, Barditch J, Gossweiler S, McNeil J, Smith P, Buldoc F, Scott R, Certa U, Broger C, Tully T: The staufen/pumilio pathway is involved in Drosophila long-term memory. Curr Biol. 2003, 13 (4): 286-296.

    PubMed  CAS  Article  Google Scholar 

  49. 49.

    Fedorowicz GM, Fry JD, Anholt RRH, Mackay TFC: Epistatic interactions between smell-impaired loci in Drosophila melanogaster. Genetics. 1998, 148 (4): 1885-1891.

    PubMed  CAS  PubMed Central  Google Scholar 

  50. 50.

    Lilly M, Kreber R, Ganetzky B, Carlson JR: Evidence That the Drosophila Olfactory Mutant Smellblind Defines a Novel Class of Sodium Channel Mutation. Genetics. 1994, 136: 1087-1096.

    PubMed  CAS  PubMed Central  Google Scholar 

  51. 51.

    Martin C, Provost E, Roux M, Bruchou C, Crauser D, Clement JL, Le Conte Y: Resistance of the honey bee, Apis mellifera to the acarian parasite Varroa destructor: behavioural and electroantennographic data. Physiol Entomol. 2001, 26 (4): 362-370.

    Article  Google Scholar 

  52. 52.

    Gramacho KP, Spivak M: Differences in olfactory sensitivity and behavioral responses among honey bees bred for hygienic behavior. Behavioral Ecology and Sociobiology. 2003, 54 (5): 472-479.

    Article  Google Scholar 

  53. 53.

    Masterman R, Ross R, Mesce K, Spivak M: Olfactory and behavioral response thresholds to odors of diseased blood differ between hygienic and non-hygienic honey bees (Apis mellifera L.). J Comp Physiol A. 2001, 187 (6): 441-452.

    PubMed  CAS  Article  Google Scholar 

  54. 54.

    Masterman R, Smith BH, Spivak M: Brood odor discrimination abilities in hygienic honey bees (Apis mellifera L.) using proboscis extension reflex conditioning. Journal of Insect Behavior. 2000, 13 (1): 87-101.

    Article  Google Scholar 

  55. 55.

    Willoughby L, Batterham P, Daborn PJ: Piperonyl butoxide induces the expression of cytochrome P450 and glutathione S-transferase genes in Drosophila melanogaster. Pest Management Science. 2007, 63 (8): 803-808.

    PubMed  CAS  Article  Google Scholar 

  56. 56.

    Pittendrigh B, Reenan R, ffrench-Constant RH, Ganetzky B: Point mutations in the Drosophila sodium channel gene para associated with resistance to DDT and pyrethroid insecticides. Mol Gen Genet. 1997, 256 (6): 602-610.

    PubMed  CAS  Article  Google Scholar 

  57. 57.

    Vogel C, Teichmann SA, Chothia C: The immunoglobulin superfamily in Drosophila melanogaster and Caenorhabditis elegans and the evolution of complexity. Development. 2003, 130 (25): 6317-6328.

    PubMed  CAS  Article  Google Scholar 

  58. 58.

    Yao B, Rakhade SN, Li QF, Ahmed S, Krauss R, Draghici S, Loeb JA: Accuracy of cDNA microarray methods to detect small gene expression changes induced by neuregulin on breast epithelial cells. Bmc Bioinformatics. 2004, 5: Art. No. 99.

    Google Scholar 

  59. 59.

    Cash AC, Whitfield CW, Ismail N, Robinson GE: Behavior and the limits of genomic plasticity: power and replicability in microarray analysis of honeybee brains. Genes Brain Behav. 2005, 4 (4): 267-271.

    PubMed  CAS  Article  Google Scholar 

  60. 60.

    Whitfield CW, Ben-Shahar Y, Brillet C, Leoncini I, Crauser D, LeConte Y, Rodriguez-Zas S, Robinson GE: Genomic dissection of behavioral maturation in the honey bee. PNAS. 2006, 103 (44): 16068-16075.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  61. 61.

    Whitfield CW, Cziko AM, Robinson GE: Gene expression profiles in the brain behavior in individual honey bees. Science. 2003, 302 (5653): 296-299.

    PubMed  CAS  Article  Google Scholar 

  62. 62.

    Ellis JDJ: The future of Varroa control: integrating current treatments with the latest advancements. Am Bee J. 2001, 141 (2): 127-131.

    Google Scholar 

  63. 63.

    Trouiller J, Milani N: Stimulation of Varroa jacobsoni Oud. oviposition with semiochemicals from honeybee brood. Apidologie. 1999, 30 (1): 3-12.

    Article  Google Scholar 

  64. 64.

    Nazzi F, Milani N, Della Vedova GD: (Z)-8-heptadecene from infested cells reduces the reproduction of Varroa destructor under laboratory conditions. J Chem Ecol. 2002, 28 (11): 2181-2190.

    PubMed  CAS  Article  Google Scholar 

  65. 65.

    Whinston ML: Development and nutrition. The biology of the honeybee. 1987, Boston, MA: Harvard University Press, 46-71.

    Google Scholar 

  66. 66.

    Winston ML: The biology of the honey bee. 1987, Cambridge, MA: Harvard University Press

    Google Scholar 

  67. 67.

    Franck P, Coussy H, Le Conte Y, Solignac M, Garnery L, Cornuet J-M: Microsatellite analysis of sperm mixture in honeybee. Insect Mol Biol. 1999, 8: 419-421.

    PubMed  CAS  Article  Google Scholar 

  68. 68.

    R Development Core Team: R: A language and environment for statistical computing. 2005, Vienna: R Foundation for Statistical Computing, 2.2.0

    Google Scholar 

  69. 69.

    Cui X, Churchill G: Statistical tests for differential expression in dCNA microarray experiments. Genome Biology. 2003, 4: 210-

    PubMed  PubMed Central  Article  Google Scholar 

  70. 70.

    Cui XQ, Hwang JTG, Qiu J, Blades NJ, Churchill GA: Improved statistical tests for differential gene expression by shrinking variance components. Biostatistics. 2005, 6 (1): 59-75.

    PubMed  Article  Google Scholar 

  71. 71.

    Whitfield CW, Band MR, Bonaldo MF, Kumar CG, Liu L, Pardinas JR, Robertson HM, Bento Soares M, Robinson GE: Annotated expressed sequence tags and cDNA microarrays for studies of Brain and Behavior in the Honey bee. Genome Res. 2002, 12: 555-566.

    PubMed  PubMed Central  Article  Google Scholar 

  72. 72.

    BEEBASE. The Honey Bee Genome Database. []

Download references


We thank C. Whitfield for helpful discussion in a previous version of the experimental design, advice in hybridization protocols and earlier statistical analysis. We also thank A. Rieux for his input in the real-time PCR experiments and J-M Bécard and T. Newman for technical assistance. We thank B. Blondin and P. Delobel for help in array's scanning and for providing access to laboratory facilities at Montpellier SupAgro. This work was funded by a European grant for beekeeping, FEOGA-convention 07-17/2004-2007.

Author information



Corresponding author

Correspondence to M Navajas.

Additional information

Authors' contributions

MN participated in conceiving the project and in designing the experimental setup, oversaw the research and primarily wrote the paper. AM participated in the planning of all experiments and performed statistical analysis. CA developed behavioral interpretation of results and actively contributed to the writing of the paper. MLM–M designed statistical analysis. GER provided support to apply honeybee microarray resources and contributed to discussion of results. JDE performed analysis of results related to immunity and contributed to discussion of overall results. SC–A performed technical work and contributed to data acquisition. DC performed technical beekeeping experiments. YLC participated in the conception and design of the study, conceived and carried out all biological experiences and participated in discussion of overall results. All authors read, edited, and approved the final manuscript.

A Migeon, C Alaux contributed equally to this work.

Electronic supplementary material

List of

Additional file 1: Apis mellifera genes with significantly different expression profiles in Varroa-parasitized and non-parasitized bees (A); in two types of bee genotypes (either sensitive or tolerant to Varroa) (B). The EST ID, the official honey bee gene ID, the gene name of the best symbol of Drosophila melanogaster (based on Flybase information) and the fold change in gene expression, are indicated for both up-regulated and down-regulated transcripts. Unassigned genes are indicated by *. (DOC 192 KB)

Additional file 2: Enrichment analysis on genes that show significant differences in bees parasitized by Varroa (A) and between two different bee genotypes (tolerant and sensitive to Varroa) (B). GO ID and GO term refer respectively to Gene Ontology ID and the associated term. (PDF 51 KB)

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and Permissions

About this article

Cite this article

Navajas, M., Migeon, A., Alaux, C. et al. Differential gene expression of the honey bee Apis mellifera associated with Varroa destructor infection. BMC Genomics 9, 301 (2008).

Download citation


  • Hygienic Behavior
  • Varroa Destructor
  • Deformed Wing Virus
  • Varroa Mite
  • Queen Mandibular Pheromone