<title>Background</title> <p>The honey bee (<italic>Apis mellifera</italic>, Insecta: Hymenoptera) has become an important model for genetic study, especially as its genome has been sequenced [##REF##17073008##1##]. It is also an important economic insect as it is the world's principal crop pollinator and honey producer [##UREF##0##2##]. These activities have been threatened by the spread of <italic>Varroa destructor </italic>(Acari: Parasitiformes) a parasite of honey bees that causes devastating harm in many countries [##UREF##1##3##]. <italic>Varroa </italic>mites are ectoparasites of honey bees, parasitizing immature and adult bees and reproducing in cells in the honeycomb that contain brood [##REF##10761588##4##]. <italic>Varroa </italic>mites or virus associated to mites impair the honey bee immune system [##UREF##2##5##] and in some cases boost the amplification of bee viruses [##REF##15897457##6##]. The mechanisms underlying the mite's suppression of bee immunity and its impacts on pathogen virulence have not been elucidated.</p> <p>There are genetic differences in the ability of honey bees to tolerate <italic>Varroa </italic>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 <italic>Varroa </italic>infection and survive despite the presence of the parasite in the hive [##UREF##3##7##, ####UREF##4##8##, ##UREF##5##9##, ##UREF##6##10####6##10##]. These differences have been attributed to a variety of factors, including grooming and hygienic behavior [##REF##10761588##4##], and differences in the timing of larval and pupal development that impact mite reproduction [##UREF##7##11##].</p> <p>Colonies of the Asian honey bee <italic>Apis cerana </italic>(the original host of <italic>V. destructor</italic>[##REF##11108385##12##]) suffer less damage from this parasite than <italic>A. mellifera </italic>in spite of the presence of the mite in the hives, and these two factors have been implicated [##UREF##8##13##]. The mechanisms underlying genetic differences for honey bee tolerance to <italic>Varroa </italic>mites are unknown. Insights into these mechanisms may lead to new molecular tools for both <italic>Varroa </italic>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 [##REF##16790040##14##, ####REF##16164603##15##, ##REF##16061818##16##, ##REF##14749722##17##, ##REF##16277749##18####16277749##18##].</p>
<title>Methods</title> <title>Honey bee colonies and sample collection</title> <p>Both, sensitive and tolerant honey bee colonies here studied, belong to non related local strains of the same <italic>A. mellifera </italic>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 <italic>Varroa </italic>despite the widespread presence of this parasite in the locale, and displayed a very low rate of parasitism [##UREF##6##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 et.al. [##UREF##15##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.</p> <p>The reproduction of <italic>Varroa </italic>mites can be affected by traits in developing bees [##UREF##16##63##,##REF##12523561##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 [##UREF##17##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 <italic>Varroa </italic>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.</p> <p>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 [##UREF##18##66##]) were collected in equal numbers from each of the study colonies.</p> <title>Sibling estimation by microsatellite genotyping</title> <p>To identify full sisters, one hundred bees per colony were genotyped (50 parasitized and 50 not parasitized) using <italic>Ap53 </italic>and <italic>A107 </italic>microsatellite loci; these loci have been shown to be effective for this purpose [##REF##10469260##67##]. DNA was extracted from 2 legs/pupa, and DNA extraction and PCR amplification were as described [##REF##10469260##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 <italic>Varroa</italic>-parasitized and non parasitized bees were analyzed in microarray experiments. This procedure minimizes the effects of intra-colonial genetic variation on gene expression.</p> <title>Microarray experimental design</title> <p>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. ##FIG##2##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.</p> <title>RNA isolation</title> <p>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.</p> <title>cDNA synthesis and array hybridizations</title> <p>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 [##REF##14551438##61##]. The slides were hybridized at 42°C and air dried by centrifugation for 2 min at 800 rpm at room temperature.</p> <title>Microarray data acquisition and statistical analysis</title> <p>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 [##REF##14551438##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.</p> <p>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 &lt; 0.05 was used to denote statistical significance).</p> <p>Statistical analyses were conducted using R version 2.2.0 [##REF##14551438##61##] and the R/maanova software package version 0.98.8 [##REF##14551438##61##,##UREF##19##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 <italic>Varroa</italic>, 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 [##REF##12702200##69##,##REF##15618528##70##]. To control the number of false-positives, the significance (nominal <italic>P</italic>-value) of genes was computed using a permutation test (pvalperm with 1000 permutations) [##REF##15618528##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 <italic>P</italic>-value is lower than 0.05. The adjusted <italic>P</italic>-value used here made it possible to control the Family Wise Error Rate (FWER)</p> <title>Functional analysis of gene expression</title> <p>Microarray ESTs were annotated as described in [##REF##17065327##60##,##REF##11932240##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 [##UREF##20##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 <italic>Drosophila melanogaster</italic>.</p> <p>We used Gene Ontology (GO) analysis to identify the biological function of the differentially expressed genes, as described in [##REF##17065327##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 [##REF##17065327##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.</p> <title>Transcript quantification by RT-qPCR</title> <p>The results of the microarray study were confirmed by measuring the expression of 4 differentially regulated genes (Table ##TAB##1##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 dH<sub>2</sub>O and 2 μl cDNA template. Three reactions were performed and means calculated for each locus and treatment. To standardize the results, a housekeeping gene (<italic>Rp49) </italic>that did not vary in expression was used as control. The forward/reverse primer sequences are indicated in Table ##TAB##1##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 [##REF##17065327##60##] method. Efficiency of the amplicons obtained for each locus was adequately high and at least 97% (Table ##TAB##1##2##).</p>
<title>Results</title> <title>Differences in bee gene expression in response to Varroa infestation and bee genotype</title> <p>Comparisons of <italic>Varroa</italic>-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 <italic>Varroa</italic>. 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 <italic>Varroa</italic>. BLAST searches indicated that this EST matches the honey bee virus (AY292384), deformed wing virus [##REF##16641291##19##].</p> <p>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 <italic>Varroa </italic>infection and bee genotype. <italic>Dlic </italic>2 is down-regulated in <italic>Varroa</italic>-parasitized bees and up-regulated in tolerant bees, whereas <italic>Strn-Mlck </italic>is down-regulated both in <italic>Varroa</italic>-parasitized and tolerant bees. The full list of genes with significant differences in expression is presented in Additional file ##SUPPL##0##1##.</p> <title>Identification of <italic>A. mellifera</italic>-regulated transcripts and assignment of biological functions</title> <p>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 <italic>Varroa </italic>and 68 associated with bee genotype had information in Gene Ontology.</p> <p>To unravel functional differences related to <italic>Varroa </italic>parasitism and bee genotype, we analyzed sets of genes for overrepresentation in functional GO categories according to biological process and molecular function (Additional file ##SUPPL##1##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 <italic>Varroa </italic>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 ##TAB##0##1##. Most of the genes affected by <italic>Varroa </italic>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.</p> <title>RT-qPCR analysis</title> <p>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 ##TAB##1##2##). Differently expressed levels of cDNAs were quantified and are shown in Fig. ##FIG##0##1##. In accord with microarray analysis, cDNAs from <italic>Varroa</italic>-parasitized bees were up regulated for genes <italic>baz </italic>and CG9520 (unassigned) (Fig. ##FIG##0##1a## and ##FIG##0##1b##); or from tolerant bees for genes <italic>Alh </italic>and <italic>Hr78 </italic>(Fig. ##FIG##0##1c## and ##FIG##0##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. ##FIG##0##1##).</p>
<title>Discussion</title> <p>Previous gene expression studies concerning honey bee immunity have mainly investigated responses to microbial pathogens [##REF##16737710##20##]. Yet physiological responses to macroparasites such as <italic>Varroa </italic>are likely to be very different from microorganisms, as has been shown recently in a study on <italic>Drosophila </italic>revealing the ability of this species to activate a systemic immune response adapted to the invader [##REF##16061818##16##]. The results described here pinpoint several genes regulated by <italic>Varroa </italic>parasitism that can be linked to honey bee responses to the presence of <italic>Varroa </italic>or to differences in bee tolerance.</p> <title>Honey bee biological response to <italic>Varroa </italic>parasitism</title> <p>The pathways suggested by this study to be associated with honey bee response to <italic>Varroa </italic>parasitism are schematized hypothetically in Fig. ##FIG##1##2##. Responses to <italic>Varroa </italic>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 <italic>Varroa </italic>parasitism, and a second category of genes that are related to cognitive impairment.</p> <title><italic>Varroa </italic>parasitism and deformed wing virus</title> <p>One of the consequences of <italic>Varroa </italic>parasitism is a decline in immune capacity which appears to induce the proliferation of viruses such as deformed wings virus in bees [##REF##15897457##6##]. The down-regulation of the autophagic-specific gene 18 (<italic>Atg18</italic>) detected here is noticeable. Autophagy is an important mechanism for innate immunity against bacteria and viruses [##REF##16099218##21##,##REF##15680321##22##]. The candidate innate immunity gene <italic>poly U binding factor 68 Kd </italic>(pUf68) [##UREF##9##23##] is also down-regulated in <italic>Varroa </italic>parasitized bees. By decreasing autophagy and immunity processes in bees, <italic>Varroa </italic>might favor the proliferation of DWV. Interestingly, the <italic>Varroa</italic>-parasitized bees displayed high levels of DWV viral RNA (about 20-fold; Additional file ##SUPPL##0##1##). The boost of DWV multiplication might cause cellular and molecular damage, and thus the observed production of genes for protein repair (<italic>Pcmt) </italic>[##REF##9061928##24##], and the labeling of proteins for degradation (<italic>Nedd8</italic>) [##REF##12231629##25##]. In contrast to these findings, we did not see a decrease on transcript abundance of immune pathway members [##REF##16737710##20##] (Evans et al., 2006) found on this array. In fact, transcripts for the gene <italic>Rab7</italic>, a plausible regulator of immunity, were up-regulated in <italic>Varroa</italic>-parasitized bees.</p> <p>The most notable symptoms of <italic>Varroa</italic>-parasitized bees are disfigured, small adults with deformed legs and wings [##REF##10761588##4##,##REF##15897457##6##]. Our results show that at the transcriptome level the presence of <italic>Varroa </italic>down-regulates two genes involved in developmental processes, <italic>slg </italic>and <italic>dlg1</italic>. It has been shown that the <italic>Drosophila sugarless </italic>(<italic>sgl</italic>) gene regulates <italic>wingless </italic>signaling [##REF##9342049##26##], which has a critical role in developmental processes. The gene <italic>dlg1 </italic>has been implicated in the control of proliferation of <italic>Drosophila </italic>imaginal discs programmed to produce adult structures at metamorphosis [##REF##8830775##27##]. Although a down-regulation of these genes could have been induced by the presence of <italic>Varroa </italic>(via the DWV virus), a specific link with the development of deformed adults of virus-infected bees would need to be further investigated.</p> <title>Cognitive impairment in bees parasitized by <italic>Varroa</italic></title> <p><italic>Varroa </italic>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 [##REF##17123087##28##], prolonged absences from the nest and a lower rate of return to the colony [##UREF##10##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 <italic>Varroa </italic>remain unknown, our results show that, the gene <italic>pale </italic>encoding tyramine hydroxylase is down-regulated in pupae parasitized by <italic>Varroa</italic>. Interestingly, the gene <italic>pale </italic>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 <italic>Dlic2 </italic>and <italic>Atg18 </italic>genes, both down-regulated in <italic>Varroa</italic>-parasitized bees, are enhancers of the <italic>blue cheese </italic>gene (<italic>bchs</italic>)[##REF##17435236##30##], which is up-regulated in tolerant bees (see below). The <italic>bchs </italic>gene has been reported as preventing progressive neural degeneration in aged flies [##REF##12598614##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.</p> <title>Behavioral resistance of bees to <italic>Varroa</italic></title> <p>Grooming behavior seems to be one of the main characters involved in the tolerance of bees to <italic>Varroa </italic>[##REF##10761588##4##]. Similarly, selection for tolerance to <italic>Varroa </italic>seems to be possible by selecting bees with a high level of hygienic behavior [##UREF##11##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 <italic>Varroa</italic>-tolerant and sensitive bees.</p> <p>A disproportionately high fraction of the genes differentially expressed between tolerant and susceptible bees are involved in the development of the nervous system (Fig. ##FIG##1##2B##). A large part of these genes are down-regulated in tolerant bees compared to sensitive bees, including the genes <italic>futsch, scratch (scrt), otk, Myosin heavy-chain-like (Mhcl), groucho </italic>(<italic>gro), kekkon</italic>-1 (<italic>kek1</italic>) [##REF##9342049##26##,##REF##15456725##33##, ####REF##12546818##34##, ##REF##15534205##35##, ##REF##8812109##36##, ##REF##8001118##37##, ##REF##7557390##38####7557390##38##]. <italic>Fringe </italic>(<italic>fng</italic>), also down regulated in <italic>Varroa</italic>-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 [##REF##17409286##39##]. Finally, a gene involved in locomotory behavior, <italic>single-minded </italic>(<italic>sim</italic>), is down-regulated in tolerant bees. <italic>Drosophila </italic>mutant <italic>sim </italic>flies are only able to walk in circles and this phenotype is due to defects in the central brain complex [##REF##12221007##40##].</p> <p>Several genes involved in neuron excitability are up-regulated in <italic>Varroa</italic>-tolerant bees compared to sensitive bees. <italic>Purity of essence </italic>(<italic>poe</italic>) (also named <italic>pushover</italic>) 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 [##REF##8846899##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) [##REF##14573707##42##]. The <italic>GluClα </italic>gene, a glutamate-gated Cl<sup>- </sup>channel specific to arthropods, which is also up-expressed in tolerant bees, is known to modulate neuronal membrane excitability [##REF##10966106##43##]. The paralytic gene (<italic>para</italic>), known to be important in the conducting of nerve action potentials in flies [##REF##6096206##44##], is also up-regulated. Seen together, these genes suggest a mechanism by which <italic>Varroa</italic>-tolerant bees are more sensitive to external stimuli than <italic>Varroa </italic>sensitive bees. Also, the Dynein heavy chain 64C gene (<italic>Dhc64c</italic>), up-regulated in tolerant bees, is required for proliferation of mushroom-body neuroblasts [##REF##11056531##45##]. Mushroom bodies have an important role in insect cognition and are known to be involved in learning and memory, particularly for smell.</p> <p>Interestingly, several genes involved in olfaction (smell impaired 21F <italic>smi21F</italic>, suppressor of the white-apricot gene <italic>su(w</italic><sup><italic>a</italic></sup><italic>)</italic>, <italic>poe</italic>, <italic>para, Rogdi</italic>) are up-regulated in <italic>Varroa</italic>-tolerant bees compared to sensitive bees [##REF##11238254##46##, ####REF##8722782##47##, ##REF##12593794##48##, ##REF##9560402##49##, ##REF##8005416##50####8005416##50##]. However, the Down syndrome cell adhesion molecule gene (<italic>Dscam</italic>) and <italic>four wheel drive </italic>(<italic>fwd</italic>) also have a role in olfaction, and they are down-regulated in <italic>Varroa</italic>-tolerant bees [##REF##12546818##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 <italic>Varroa </italic>infested cells. Observations made on two earlier generations of the bee stocks here studied have been shown that <italic>Varroa</italic>-tolerant bees have better ability to detect the mite [##UREF##12##51##]. Previous studies showed that hygienic bees have a higher olfactory sensitivity and responsiveness than non-hygienic bees [##UREF##13##52##, ####REF##11548991##53##, ##UREF##14##54####14##54##]. 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 <italic>Varroa</italic>. If these gene candidate pathways for <italic>Varroa </italic>behavioral tolerance can be further confirmed, the results suggest that they already set up during the pupae stage.</p> <p>There were differences in expression between tolerant and sensitive bees for genes involved in increased resistance to toxins like Ahcy13 (<italic>Adenosylhomocysteinase at 13</italic>) involved in detoxification [##REF##17514638##55##], and <italic>para </italic>involved in resistance to insecticide pyrethroids and DDT (Dichloro-Diphenyl-Trichloroethane) [##REF##9435785##56##]. Two genes linked to immunity, <italic>Dsam </italic>and <italic>otk</italic>, members of the Immunoglobulin gene superfamily in <italic>Drosophila </italic>[##REF##14623821##57##] were down-regulated in <italic>Varroa</italic>-tolerant bees, while one galectin-family gene (an apparent ortholog to Drosophila CG32226) appears to be upregulated.</p> <title>Low fold differences in gene expression</title> <p>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 [##REF##15272935##58##]; 87% of a set of parasite-specific genes displayed changes 2-fold or less in <italic>D. melanogaster </italic>challenged to a protozoan parasite [##REF##14749722##17##]. Similarly, about 60% of genes of the lepidopteran, <italic>Spodoptera frugiperda</italic>, revealed around 0.5-fold changes in the transcript levels associated with virus infection [##REF##16790040##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 [##REF##16277749##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 <italic>Varroa </italic>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 [##REF##14573707##42##,##REF##15924558##59##, ####REF##17065327##60##, ##REF##14551438##61####14551438##61##].</p>
<title>Conclusion</title> <p>This work is the first step towards understanding the genomic responses of honey bees to <italic>Varroa </italic>parasitism. It demonstrates that honey bee pupae exhibit differences in gene expression associated either with the presence of <italic>Varroa</italic>, or with tolerance to this parasite. These results highlight the potential importance of behavioral mechanisms of response to <italic>Varroa </italic>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.</p>
<p>This is an Open Access article distributed under the terms of the Creative Commons Attribution License (<ext-link ext-link-type="uri" xlink:href="http://creativecommons.org/licenses/by/2.0"/>), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.</p> <title>Background</title> <p>The parasitic mite, <italic>Varroa destructor</italic>, is the most serious pest of the western honey bee, <italic>Apis mellifera</italic>, 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 <italic>Varroa </italic>infestation induces changes in <italic>Apis mellifera </italic>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 <italic>Varroa </italic>parasitism.</p> <title>Results</title> <p>We explored the transcriptional response to mite parasitism in two genetic stocks of <italic>A. mellifera </italic>which differ in susceptibility to <italic>Varroa</italic>, 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 <italic>Varroa</italic>-susceptible and two colonies of <italic>Varroa</italic>-tolerant bees. We identified a set of 148 genes with significantly different patterns of expression: 32 varied with the presence of <italic>Varroa</italic>, 116 varied with bee genotype, and 2 with both. <italic>Varroa </italic>parasitism caused changes in the expression of genes related to embryonic development, cell metabolism and immunity. Bees tolerant to <italic>Varroa </italic>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 <italic>Varroa</italic>; differences in olfaction may be related to increased grooming and hygienic behavior, important behaviors known to be involved in <italic>Varroa </italic>tolerance.</p> <title>Conclusion</title> <p>These results suggest that differences in behavior, rather than in the immune system, underlie <italic>Varroa </italic>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.</p>
<title>Authors' contributions</title> <p>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.</p> <title>Supplementary Material</title>
<title>Acknowledgements</title> <p>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.</p>
<fig position="float" id="F1"><label>Figure 1</label><caption><p><bold>cDNA levels of four transcripts shown by microarray analysis to be significantly regulated in <italic>Apis mellifera</italic></bold>. In a and b, samples of bees infested with <italic>Varroa </italic>(hatched bars) are compared with bees free of <italic>Varroa </italic>(white bars). In c and d: <italic>Varroa</italic>-sensitive bees (hatched bars) are compared with <italic>Varroa</italic>-tolerant bees (white bars). All the values shown are mean ± SE.</p></caption></fig> <fig position="float" id="F2"><label>Figure 2</label><caption><p><bold>Hypothetical pathways and models of honey bee responses to <italic>Varroa</italic>-parasitism (A) and the bee tolerant genotype (B)</bold>. Arrows and dashes indicate positive and negative regulation, respectively. Gene names in bold are up-regulated. In A) one of the consequences of the <italic>Varroa </italic>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 (<italic>Pcmt</italic>) and the labelling of proteins for degradation (<italic>Nedd8</italic>). In addition, regulated genes that might be affected by the presence of the DWV are indicated. Mites might decrease the production of dopamine (<italic>ple</italic>) and inhibit genes known for indirectly preventing neural degeneration in aged adults, which could explain the cognitive impairment often observed in adults parasitized by <italic>Varroa</italic>. In B) different genes can be associated to behavioral tolerance to <italic>Varroa</italic>. 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 <italic>Varroa</italic>-tolerance. See text for a full discussion on the genes involved in the pathways presented.</p></caption></fig> <fig position="float" id="F3"><label>Figure 3</label><caption><p><bold>Microarray experimental design</bold>. For all experiments, arrows indicate microarray hybridizations (arrow tail, Cy3-labeled sample; arrowhead, Cy5-labeled sample). <italic>Varroa </italic>parasitized (+) and non-parasitized (-) full-sister pupae, from two different genetic stocks: one susceptible (S) and one tolerant (T) to <italic>Varroa </italic>were compared. Dye swaps were made for all comparisons.</p></caption></fig>
<table-wrap position="float" id="T1"><label>Table 1</label><caption><p>Functional clustering based on Gene Ontology</p></caption><table frame="hsides" rules="groups"><thead><tr><td align="left" colspan="4"><bold>(A) Varroa</bold></td></tr></thead><tbody><tr><td/><td align="center"><italic>P </italic>value</td><td align="center">Up-regulated</td><td align="center">Down-regulated</td></tr><tr><td colspan="4"><hr/></td></tr><tr><td align="left"><bold>Biological process </bold>(unknown)</td><td/><td align="center">6 (9)</td><td align="center">10 (6)</td></tr><tr><td align="left">cell organization and biogenesis</td><td align="center">0.000229</td><td align="center">2</td><td align="center">3</td></tr><tr><td align="left">protein metabolism</td><td align="center">0.004808</td><td align="center">4</td><td align="center">5</td></tr><tr><td align="left">cellular metabolism</td><td align="center">0.008242</td><td align="center">4</td><td align="center">6</td></tr><tr><td align="left">protein modification</td><td align="center">0.035714</td><td align="center">2</td><td align="center">1</td></tr><tr><td align="left"><bold>Molecular function </bold>(unknown)</td><td/><td align="center">5 (9)</td><td align="center">9 (8)</td></tr><tr><td align="left">transferase activity</td><td align="center">0.021978</td><td align="center">3</td><td align="center">0</td></tr><tr><td align="left">catalytic activity</td><td align="center">0.041958</td><td align="center">3</td><td align="center">3</td></tr><tr><td colspan="4"><hr/></td></tr><tr><td align="left" colspan="4"><bold>(B) Bee phenotype</bold></td></tr><tr><td colspan="4"><hr/></td></tr><tr><td/><td align="center"><italic>P </italic>value</td><td align="center">Up-regulated</td><td align="center">Down-regulated</td></tr><tr><td colspan="4"><hr/></td></tr><tr><td align="left"><bold>Biological process </bold>(unknown)</td><td/><td align="center">24 (20)</td><td align="center">24 (31)</td></tr><tr><td align="left">regulation of physiological process</td><td align="center">3.279e-08</td><td align="center">3</td><td align="center">6</td></tr><tr><td align="left">regulation of transcription</td><td align="center">2.650e-09</td><td align="center">3</td><td align="center">5</td></tr><tr><td align="left">Axonogenesis</td><td align="center">5.139e-06</td><td align="center">0</td><td align="center">4</td></tr><tr><td align="left">cell division</td><td align="center">0.000231</td><td align="center">0</td><td align="center">3</td></tr><tr><td align="left">Localization</td><td align="center">1.450e-05</td><td align="center">7</td><td align="center">2</td></tr><tr><td align="left">axon guidance</td><td align="center">5.782e-05</td><td align="center">0</td><td align="center">3</td></tr><tr><td align="left">Phosphorylation</td><td align="center">7.709e-05</td><td align="center">1</td><td align="center">3</td></tr><tr><td align="left">cell organization and biogenesis</td><td align="center">0.000155</td><td align="center">5</td><td align="center">2</td></tr><tr><td align="left">nervous system development</td><td align="center">7.359e-05</td><td align="center">1</td><td align="center">4</td></tr><tr><td align="left">protein modification</td><td align="center">0.000360</td><td align="center">3</td><td align="center">1</td></tr><tr><td align="left"><bold>Molecular function </bold>(unknown)</td><td/><td align="center">20 (24)</td><td align="center">26 (29)</td></tr><tr><td align="left">ATPase activity</td><td align="center">5.428e-05</td><td align="center">2</td><td align="center">1</td></tr><tr><td align="left">nucleic acid binding</td><td align="center">3.480e-07</td><td align="center">4</td><td align="center">5</td></tr><tr><td align="left">methyltransferase activity</td><td align="center">5.428e-05</td><td align="center">1</td><td align="center">2</td></tr><tr><td align="left">pyrophosphatase activity</td><td align="center">4.720e-06</td><td align="center">3</td><td align="center">1</td></tr><tr><td align="left">ion transporter activity</td><td align="center">5.428e-05</td><td align="center">2</td><td align="center">1</td></tr><tr><td align="left">ubiquitin-protein ligase activity</td><td align="center">5.428e-05</td><td align="center">2</td><td align="center">1</td></tr><tr><td align="left">kinase activity</td><td align="center">0.001086</td><td align="center">2</td><td align="center">1</td></tr><tr><td align="left">catalytic activity</td><td align="center">2.904e-06</td><td align="center">6</td><td align="center">8</td></tr><tr><td align="left">transmembrane receptor protein kinase activity</td><td align="center">5.428e-05</td><td align="center">0</td><td align="center">3</td></tr><tr><td align="left">cytoskeletal protein binding</td><td align="center">0.000543</td><td align="center">0</td><td align="center">3</td></tr><tr><td align="left">transferase activity</td><td align="center">1.339e-09</td><td align="center">3</td><td align="center">7</td></tr><tr><td align="left">ion binding</td><td align="center">1.101e-05</td><td align="center">2</td><td align="center">2</td></tr><tr><td align="left">signal transducer activity</td><td align="center">0.002336</td><td align="center">0</td><td align="center">4</td></tr></tbody></table></table-wrap> <table-wrap position="float" id="T2"><label>Table 2</label><caption><p>Characteristics of the RT-qPCR analysis of differentially expressed genes of <italic>Apis mellifera</italic></p></caption><table frame="hsides" rules="groups"><thead><tr><td align="left">Product name</td><td align="center">Honey bee ID</td><td align="left">Forward primer sequences</td><td align="left">Reverse primer sequence</td><td align="center">Amplicon size (bp)</td><td align="right">Amplification efficiency<sup>a</sup>(R<sup>b</sup>)</td></tr></thead><tbody><tr><td align="left"><italic>Krh1</italic></td><td align="center">GB14867</td><td align="left">ACTCATCAGTTGTTGGTTCTCCTC</td><td align="left">TCGTTTGGCTCTTCAGTCTTGTG</td><td align="center">118</td><td align="right">2.11 (0.99)</td></tr><tr><td align="left">non determined</td><td align="center">GB18554</td><td align="left">TCACACCGATATTCTCATCAAAGG</td><td align="left">CTTGTCATTCTTGTTCTCCGATTG</td><td align="center">112</td><td align="right">1.94 (0.98)</td></tr><tr><td align="left"><italic>Hr78</italic></td><td align="center">GB18358</td><td align="left">TGACGAAGTTTAGTTGCTGCTATG</td><td align="left">TGTTGTTCCCTATGATCTCTGTCC</td><td align="center">107</td><td align="right">1.98 (0.99)</td></tr><tr><td align="left"><italic>Alh</italic></td><td align="center">GB17400</td><td align="left">ACTTGTGGTAATGCTGGCTGAC</td><td align="left">AACGAACGAAGGAAAGGAATAACG</td><td align="center">129</td><td align="right">2.08 (0.96)</td></tr><tr><td align="left"><italic>baz</italic></td><td align="center">GB10346</td><td align="left">ACCAGGAACAAGCGAGTCAGAAG</td><td align="left">ACCAGGAACAAGCGAGTCAGAAG</td><td align="center">112</td><td align="right">2.02 (0.99)</td></tr></tbody></table></table-wrap>
<supplementary-material content-type="local-data" id="S1"><caption><title>Additional file 1</title><p><bold>List of <italic>Apis mellifera </italic>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)</bold>. The EST ID, the official honey bee gene ID, the gene name of the best symbol of <italic>Drosophila melanogaster </italic>(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 *.</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="S2"><caption><title>Additional file 2</title><p><bold>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)</bold>. GO ID and GO term refer respectively to Gene Ontology ID and the associated term.</p></caption></supplementary-material>
<table-wrap-foot><p>Functional clustering based on Gene Ontology of genes that show significant differences are indicated for bees parasitized by Varroa (A), and for Varroa-tolerant bees (B). Gene Ontology terms are non-exclusive. Expression patterns within these sets of genes are represented (number of genes that are significantly up- or down-regulated). Genes with unknown biological process or molecular function are in brackets.</p></table-wrap-foot> <table-wrap-foot><p>Honey bee IDs refer to the cDNAs selected from the EST microarray experiment (see text).</p><p><italic>Krh1 </italic>is the housekeeping gene. Forward and reverse primer sequences and the PCR product length are indicated.</p><p><sup>a</sup>Amplification efficiencies were calculated from the slope of standard curves as E = 10 [-1/slope] – 1.100% efficiency corresponds to an amplification efficiency of 1; <sup>b</sup>Regression coefficient of linear standard curve.</p></table-wrap-foot>
<graphic xlink:href="1471-2164-9-301-1"/> <graphic xlink:href="1471-2164-9-301-2"/> <graphic xlink:href="1471-2164-9-301-3"/>
<media xlink:href="1471-2164-9-301-S1.doc" mimetype="application" mime-subtype="msword"><caption><p>Click here for file</p></caption></media> <media xlink:href="1471-2164-9-301-S2.pdf" mimetype="application" mime-subtype="pdf"><caption><p>Click here for file</p></caption></media>
[{"surname": ["Burgett", "Rucker", "Thurman"], "given-names": ["M", "RR", "WN"], "article-title": ["Economics and honey bee pollination markets"], "source": ["Am Bee J"], "year": ["2004"], "volume": ["144"], "fpage": ["269"], "lpage": ["271"]}, {"surname": ["Griffiths", "Bowman"], "given-names": ["DA", "CE"], "article-title": ["World distribution of the mite "], "italic": ["Varroa jacobsoni"], "source": ["Bee World"], "year": ["1981"], "volume": ["62"], "fpage": ["154"], "lpage": ["163"]}, {"surname": ["Gregory", "Evans", "Rinderer", "De Gruzman"], "given-names": ["PG", "JD", "TE", "L"], "article-title": ["Conditional immune-gene suppression of honeybees parasitized by "], "italic": ["Varroa "], "source": ["J Insect Scien"], "year": ["2005"], "volume": ["5"], "fpage": ["1"], "lpage": ["5"]}, {"surname": ["Fries", "Camazine", "Sneyd"], "given-names": ["I", "S", "J"], "article-title": ["Population-Dynamics of "], "italic": ["Varroa jacobsoni "], "source": ["Bee World"], "year": ["1994"], "volume": ["75"], "fpage": ["5"], "lpage": ["28"]}, {"surname": ["Seeley"], "given-names": ["TD"], "article-title": ["Honey bees of the Arnot Forest: a population of feral colonies persisting with "], "italic": ["Varroa destructor "], "source": ["Apidologie"], "year": ["2007"], "volume": ["38"], "fpage": ["19"], "lpage": ["29"]}, {"surname": ["Wood", "Hardin", "Lee"], "given-names": ["M", "B", "J"], "article-title": ["Varroa-tolerant bees keep hives buzzing"], "source": ["Agr Res Mag"], "year": ["1999"], "volume": ["47"], "fpage": ["10"], "lpage": ["12"]}, {"surname": ["Le Conte", "de Vaublanc", "Crauser", "Jeanne", "Rousselle", "B\u00e9card"], "given-names": ["Y", "G", "D", "F", "JC", "JM"], "article-title": ["Honey bee colonies that have survived "], "italic": ["Varroa destructor"], "source": ["Apidologie"], "year": ["2007"], "volume": ["38"]}, {"surname": ["Moritz"], "given-names": ["FA"], "article-title": ["Heritability of the postcapping stage in "], "italic": ["Apis mellifera "], "source": ["The Journal of Heredity"], "year": ["1985"], "volume": ["76"], "fpage": ["267"], "lpage": ["270"]}, {"surname": ["Peng", "Y", "S", "L"], "given-names": ["YS", "F", "X", "G"], "article-title": ["The resistance Mecanism of the Asian Honey Bee, Apis cerana Fabr., to an Ectoparasitic Mite, Varroa jacobsoni Oudemans"], "source": ["J Invert Path"], "year": ["1987"], "volume": ["49"], "fpage": ["54"], "lpage": ["60"]}, {"surname": ["Foley", "O'Farrell"], "given-names": ["E", "PH"], "article-title": ["Functional dissection of an innate immune response by a genome-wide RNAi screen"], "source": ["PLoS Biology"], "year": ["2004"], "volume": ["2"], "fpage": ["1091"], "lpage": ["1106"]}, {"surname": ["Kralj", "Fuchs"], "given-names": ["J", "S"], "article-title": ["Parasitic "], "italic": ["Varroa destructor ", "Apis mellifera "], "source": ["Apidologie"], "year": ["2006"], "volume": ["37"], "fpage": ["577"], "lpage": ["587"]}, {"surname": ["Harbo", "Harris"], "given-names": ["JR", "JW"], "article-title": ["Suppressed mite reproduction explained by the behavior of adult bees"], "source": ["J Apic Res"], "year": ["2005"], "volume": ["44"], "fpage": ["21"], "lpage": ["23"]}, {"surname": ["Martin", "Provost", "Roux", "Bruchou", "Crauser", "Clement", "Le Conte"], "given-names": ["C", "E", "M", "C", "D", "JL", "Y"], "article-title": ["Resistance of the honey bee, "], "italic": ["Apis mellifera ", "Varroa destructor"], "source": ["Physiol Entomol"], "year": ["2001"], "volume": ["26"], "fpage": ["362"], "lpage": ["370"]}, {"surname": ["Gramacho", "Spivak"], "given-names": ["KP", "M"], "article-title": ["Differences in olfactory sensitivity and behavioral responses among honey bees bred for hygienic behavior"], "source": ["Behavioral Ecology and Sociobiology"], "year": ["2003"], "volume": ["54"], "fpage": ["472"], "lpage": ["479"]}, {"surname": ["Masterman", "Smith", "Spivak"], "given-names": ["R", "BH", "M"], "article-title": ["Brood odor discrimination abilities in hygienic honey bees ("], "italic": ["Apis mellifera "], "source": ["Journal of Insect Behavior"], "year": ["2000"], "volume": ["13"], "fpage": ["87"], "lpage": ["101"]}, {"surname": ["Ellis"], "given-names": ["JDJ"], "article-title": ["The future of Varroa control: integrating current treatments with the latest advancements"], "source": ["Am Bee J"], "year": ["2001"], "volume": ["141"], "fpage": ["127"], "lpage": ["131"]}, {"surname": ["Trouiller", "Milani"], "given-names": ["J", "N"], "article-title": ["Stimulation of "], "italic": ["Varroa jacobsoni "], "source": ["Apidologie"], "year": ["1999"], "volume": ["30"], "fpage": ["3"], "lpage": ["12"]}, {"surname": ["Whinston"], "given-names": ["ML"], "article-title": ["Development and nutrition"], "source": ["The biology of the honeybee"], "year": ["1987"], "publisher-name": ["Boston, MA: Harvard University Press"], "fpage": ["46"], "lpage": ["71"]}, {"surname": ["Winston"], "given-names": ["ML"], "source": ["The biology of the honey bee"], "year": ["1987"], "publisher-name": ["Cambridge, MA: Harvard University Press"]}, {"collab": ["R Development Core Team"], "source": ["R: A language and environment for statistical computing"], "year": ["2005"], "edition": ["2.2.0"], "publisher-name": ["Vienna: R Foundation for Statistical Computing"]}, {"article-title": ["BEEBASE. The Honey Bee Genome Database"]}]
BMC Genomics. 2008 Jun 25; 9:301