Immune related gene expression in worker honey bee (Apis mellifera carnica) pupae exposed to neonicotinoid thiamethoxam and Varroa mites (Varroa destructor)

Varroa destructor is one of the most common parasites of honey bee colonies and is considered as a possible co-factor for honey bee decline. At the same time, the use of pesticides in intensive agriculture is still the most effective method of pest control. There is limited information about the effects of pesticide exposure on parasitized honey bees. Larval ingestion of certain pesticides could have effects on honey bee immune defense mechanisms, development and metabolic pathways. Europe and America face the disturbing phenomenon of the disappearance of honey bee colonies, termed Colony Collapse Disorder (CCD). One reason discussed is the possible suppression of honey bee immune system as a consequence of prolonged exposure to chemicals. In this study, the effects of the neonicotinoid thiamethoxam on honey bee, Apis mellifera carnica, pupae infested with Varroa destructor mites were analyzed at the molecular level. Varroa-infested and non-infested honey bee colonies received protein cakes with or without thiamethoxam. Nurse bees used these cakes as a feed for developing larvae. Samples of white-eyed and brown-eyed pupae were collected. Expression of 17 immune-related genes was analyzed by real-time PCR. Relative gene expression in samples exposed only to Varroa or to thiamethoxam or simultaneously to both Varroa and thiamethoxam was compared. The impact from the consumption of thiamethoxam during the larval stage on honey bee immune related gene expression in Varroa-infested white-eyed pupae was reflected as down-regulation of spaetzle, AMPs abaecin and defensin-1 and up-regulation of lysozyme-2. In brown-eyed pupae up-regulation of PPOact, spaetzle, hopscotch and basket genes was detected. Moreover, we observed a major difference in immune response to Varroa infestation between white-eyed pupae and brown-eyed pupae. The majority of tested immune-related genes were upregulated only in brown-eyed pupae, while in white-eyed pupae they were downregulated.


Introduction
Pollinators are crucial in almost all terrestrial ecosystems, especially in those dominated by agriculture. Among pollinators in Slovenia the most important is the Carniolan honey bee, Apis mellifera carnica [1]. Unfortunately, high mortality rates in overwintering colonies have been documented over the last few years in North America and Europe [2,3]. Besides the impact this has on ecological balance, such declines also have a great economic impact. Due to these massive losses, multiple candidate factors have been studied, including exposure to pesticides and pathogens [3].
The role of pesticides in honey bee losses, and their sub-lethal and synergistic impact, have been the subject of an increasing number of studies. Adult worker honey bees are repeatedly exposed to pesticides during the collection of pollen and nectar. By feeding developing honey bees contaminated food, entire colonies can be exposed to pesticides [4,5]. Such sub-lethal exposure to pesticides can influence their physiology and behavior [6] and can cause changes in immune response and detoxification mechanisms [7][8][9].
In particular, neonicotinoid pesticides are currently receiving special attention [10]. These insecticides are chemically similar to nicotine, thus they act antagonistically to insect nicotine acetylcholine receptors [11]. Consequently, neonicotinoids may reduce the learning performance of honey bees, making them more susceptible to pathogens and other factors that result in colony loss [12,13]. The neonicotinoid thiamethoxam is commonly used in agriculture and has been found as residue in honey bee colonies [4]. Recent studies demonstrated that even the lowest exposure dose of thiamethoxam elicits sub-lethal deleterious effects in adult honey bees [14] and toxicity of thiamethoxam on honey bee larvae [15]. In addition to pesticides, pathogenic infections can also contribute to colony mortality [1,7].
Varroa destructor is one of the most widespread honey bee parasites and vector of viral and other honey bee pathogens [16]. This ectoparasitic mite is one of rare adult honey bee parasite that also affects brood [17]. Several studies [9,16,[18][19][20][21][22] show that Varroa mites influence immune response, but results are contradictory. It has been suggested that Varroa parasitism suppresses honey bee immunity leaving them vulnerable to a variety of viruses [18][19][20][21][22][23]. While other studies of Varroa infestation in developing stages of worker bees showed increased expression of immune-related genes [9,16,24,25]. As findings are hard to interpret, we are still far from understanding the interactions between the parasite and the honey bee immune system.
The immune system of honey bees possesses obvious orthologues for the major members of the immune pathways consisting of: Toll (transmembrane signal transducing pathway), Imd (immune deficiency), JNK (intracellular signaling pathways) and JAK/STAT (Janus kinase/signal transducers and activators of transcription). Although there have been several studies on honey bee immune responses to pathogens in combination with pesticides [7,9,16,24], the effect of thiamethoxam in combination with Varroa infestation on the immune system has not yet been studied. The present study focuses on the effect of thiamethoxam and Varroa mites on honey bee immune-related gene expression in two developmental stages (white-eyed and brown-eyed pupae). According to previous studies describing the ontogeny of honey bee immune system [17,26] our hypothesis was that thiamethoxam consumed in larval stage will have an influence on Varroa induced immune response in pupae.

Material and methods
The honey bees used in this study were from Apis mellifera carnica colonies maintained at the Biotechnical Faculty, University of Ljubljana, Slovenia. Experiments were conducted in the Neuroethology Laboratory of the Biotechnical Faculty, University of Ljubljana. Colonies did not present visible symptoms of any known brood disease, e.g. American foulbrood, European foulbrood, Chalkbrood, or Stonebrood, which was also confirmed with RT-qPCR.

Experimental setup
Eight mixed honeybee colonies were established from two different honey bee colonies A and B. This was accomplished in two steps. First, the two colonies were mixed, and separated into two hives, original queens were added to each colony. One colony was deliberately infested with Varroa mites by the addition of two brood frames from a heavily Varroa infested colony (colony C; more than 5% of 50 randomly sampled pupae infested). Varroa mite infestation was determined by brood sampling. Honey bee workers from both original colonies, A and B, were distributed, approximately equally, to the newly established pre-experimental colonies. The aim of this step was to establish two pre-experimental colonies that were a balanced mixture of brood and adults from colonies A and B. Each pre-experimental colony was formed from two unsealed brood frames from original colony A and two unsealed brood frames from original colony B. In both cases two additional frames with honey was added. Two additional frames from the colony C were not used to setup colonies in the next steps. After 10 days, each brood frame from pre-experimental colonies A and B were split into a separate nucs, establishing eight starter colonies for two biological repetitions. To increase survival of the experimental colonies, an additional frame of unsealed brood from an original colony, either from A or B, was added to each experimental colony. Finally every experimental colony had one experimental frame and one "helper" frame of the brood and additional frame for feeding with prepared protein sugar patties. To ensure normal function of all colonies, an unfertilized queen of the same origin was added to each colony. Queens have successfully mated until the treatment procedure (see below) and they started laying eggs. The infested and non-infested colonies were placed 5-6 m apart inside of tide bushy area to reduce chance for reinvasion (drifting) during the experiment. The treatment procedure started immediately afterwards (day 0).
Each day the colonies were fed with a 50 g protein cake (4 g chicken egg yolk, 8 g yeast, 13 g tofu, 72 g powdered sugar in 1 ml water and 2 g olive oil; all ingredients were the products of organic farming). Four colonies, two infested with Varroa mites and two non-infested, received protein cake with thiamethoxam (10 μg/kg cake; Actara1, Syngenta). The concentration of thiamethoxam in the cake was chosen according to Pilling et al. [27] where evaluated residues of thiamethoxam in pollen collected from honey bees after foraging on flowering seed treated maize were between 1 and 7 μg/kg. Nurse bees collect small portions of protein cake and mix it with honey. This mixture is then normally stored for a short time in the nurse bee crops only for the purpose of the delivery to the larvae and it is not transported to the midgut where the digestion is taking place involving a considering amount of enzymes [28]. To the other four honey bee colonies, two infested with Varroa mites and two non-infested, cakes without thiamethoxam were provided. Freshly made cakes were replaced every second day and the consumption of the cake was recorded. On day ten, capped brood cells were opened and 12 white-eyed pupae and 12 brown-eyed stage worker pupae at the start of cuticle pigmentation were collected from each colony. According to COLOSS BEEBOOK [29] white-eyed pupae represent a developmental stage of 0-5 days after cell capping and brown-eyed pupae 8-12 days after the capping. In the colonies that were deliberately infested with Varroa mites only parasitized pupae were sampled. A pupa was considered to be parasitized if a reproductive Varroa was found inside the brood cell together with the pupa. From the non-infested colonies only the pupae without Varroa mites were used for further analysis. The samples were collected and snap frozen using N 2 and stored at -80˚C until RNA extraction. The set-up of treatments is summarized in S1 Table. RNA isolation and cDNA synthesis Individual white-or brown-eyed pupae were ground in liquid nitrogen and total RNA was extracted using TRIZOL1 (Ambion1, Life Technologies) and then purified with PureLink1 RNA Mini Kit columns (Ambion1, Life Technologies). RNA concentrations were determined using NanoVue (GE Healthcare, Waukesha, USA) and eight samples of white-or brown-eyed pupae from each biological group with best quality RNA were selected for further analysis. Residual DNA was removed by incubating 1 μg of RNA with 1 U of RNase-free DNase I (Fermentas, Germany) for 30 min at 37˚C. High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, USA) was used to synthesize complementary DNA (cDNA), according to the manufacturer's instructions.

Real-time PCR
Primer pairs used for candidate genes (related to immune response traits) were those reported in Gregorc et al. [9] and Evans [30] (Table 1). For quantitative real-time PCR, 10 μL reactions were prepared, containing 5 μL of Fast Start Universal SYBR Green Master (ROX) (Roche Diagnostics GmbH, Germany), 250 nM of forward and reverse primer, DEPC treated water and 5 ng of cDNA. Amplification of targets was performed with ViiA7 System (Applied Biosystems, USA) and analyzed with QuantStudio™ Real-Time PCR Software. For the experimental run the following cycle profile was used: first cycle at 95˚C for 10 min and 40 cycles at 95˚C for 20 s, 20 s at Tm of each primer pair and 72˚C for 20s, followed by last cycle for 15 s at 95˚C, 60s at Tm and 15 s at 95˚C for dissociation curve. Reactions for quantitative real-time PCR were carried out in 384-well plates (MicroAmp1, Life Technologies). Each experiment contained three no-template controls and test samples performed in duplicates. Gene expression was analyzed for 17 immune-related genes. As reference genes we used Rp49, RPS5 and Tbp-af according to its stable expression in our and previous studies [8,31]. Gene expression values of non-treated group were used for gene expression calibration. For each gene the level of gene expression was calculated using the method described by Pfaffl, where the relative expression ratio between treated and non-treated group is based on PCR efficiency [32]. The results demonstrating the level of gene expression are shown in log 2 scale. All collected samples were also tested for the most common honey bee-pathogen targets, including DWV, using RT-qPCR, performed as described before [8]. Gene expression was assayed for 10 honey bee-pathogen targets (S2 Table). The expression of pathogen genes of all the experimental samples were evaluated by comparing threshold cycle (CT) values between treatment groups.

Statistical analysis
To analyze the effects of Varroa mite infestation, thiamethoxam treatment and their interaction on gene expression, we used linear model for fixed effects. The estimation of least squares means followed by Dunnett's post hoc test was used for pairwise comparisons among the treatment groups. The assumption of normal distribution was tested and met via examination of the residuals (coefficients of skewness and kurtosis). Differences between levels were considered to be statistically significant if the p-value was equal to or less than 0.05. All statistical analyses and plotting were carried out using R software version 3.4.1 [33] with relevant libraries (lsmeans, moments, ggplot2) [34][35][36].
Pathway analysis of honey bee immune gene expression data PathVisio, the graphical editor for biological pathways [37], was used for visualization of honey bee immune pathways. We have used the WikiPathways plugin for PathVisio to search for immune pathways in insects in the online pathway database. As a base, we used the Drosophila melanogaster pathway, (WP3830-Toll, IMD, and JAK/STAT Pathways for Immune Response to Pathogens). We have modified selected pathway diagrams: 1) all the elements presenting genes which were not annotated in the Apis mellifera genome (version Amel_4.5) were removed, 2) JNK pathway was added, 3) effector genes, such as defensin1, defensin2, and abaecin were added. Thus we developed honey bee specific Toll, IMD, JAK/STAT and JNK pathway diagrams.
The honey bee pathway diagram was used for differentially expressed immune-related genes analysis. Using the ggplot library for R, we plotted a heatmap (S1 Fig), where colors indicate the average mRNA levels of Varroa treated groups compared to average levels of mRNA in groups simultaneously treated with thiamethoxam and infested with Varroa. We mapped the colors from the scale of the heatmap to the elements of pathway diagrams (S1 Fig).

Immune-related gene expression in white-eyed pupae
In white-eyed pupae, only four genes were upregulated in at least one treatment, while 13 immune-related genes were downregulated in all cases. Thiamethoxam significantly upregulated genes defensin-1 (3.94), abaecin (3.26) and dredd (1.38), while 13 genes were significantly downregulated. Among the most downregulated were lysozyme-2 (0.13), cactus (0.44) and dorsal-1 (0.44) (Fig 1). In Varroa infested groups, abaecin (6.30), dredd (1.24) and defensin-1 (30.53) had statistically significant upregulation while the most downregulated genes were lysozyme-2 (0.13) and dorsal-1 (0.40) (Fig 2). In groups simultaneously exposed to thiamethoxam and infested with Varroa, we can see a similar gene expression pattern as in the other two groups, but upregulation of dreed was not statistically significant. Regardless of treatment, all white-eyed pupae transmembrane proteins were downregulated, and no immune-related pathways were activated, except part of the Imd pathway (dredd). The only effector genes that were upregulated were abaecin and defensin-1 (Fig 1).

Immune-related gene expression in brown-eyed pupae
In brown-eyed pupae, thiamethoxam caused mild effects in gene expression resulting in statistically significant downregulation of only spaetzle (Fig 1). In Varroa infested brown-eyed pupae, gene expression pattern was completely changed compared to white-eyed pupae. The majority of genes tested were upregulated and only two genes, PPOact and spaetzle, were downregulated (Fig 1). The highest statistically significant upregulation was detected for dorsal-1 (2.39). In Varroa infested pupae, thiamethoxam treatment significantly changed the expression of ten genes (Fig 2). The pathway diagrams show that all transmembrane proteins are upregulated, which resulted in activation of all four immune-related pathways (Toll, JAK/ STAT, IMD and JNK) and effector genes (defensin-2, lysozyme-2, and PPOact) (Fig 1).

The difference in the mean expression of each gene between groups infested with Varroa and group simultaneously infested with Varroa and treated with thiamethoxam
We compared differences in the mean expression of each gene between Varroa-infested colonies and those that were Varroa-infested and treated with thiamethoxam. In infested whiteeyed pupae we observed that consumption of thiamethoxam in larval stage led to downregulation of spaetzle, abaecin and defensin-1, but it also caused higher expression of gene lysozyme-2 ( Fig 3A). Expression of spaetzle, hopscotch, basket and PPOact genes was statistically increased in brown-eyed pupae with addition of thiamethoxam to Varroa infested group (Fig 3B). The only gene on which thiamethoxam had no effect in both white-and brown-eyed infested pupae was defensin-2.

Discussion
The strength of honey bee colony strongly depends on the brood so studies regarding the effects of various neonicotinoids on larvae should be performed [10]. It has been revealed that even sub-lethal concentrations of thiamethoxam can affect the larval development of honey bee [15]. In different studies the effect of neonicotinoids and Varroa infestation has been studied but there is limited information about their combined effect. In our study, the majority of immune-related genes were downregulated in white-eyed pupae in all treatments. In Varroainfested white-eyed pupae the only significantly upregulated genes were two antimicrobial effectors, abaecin and defensin-1 (Fig 1). Antimicrobial peptides (AMPs) are evolutionary important components of insect immune systems. They can be activated and delivered in a short period of time to the site of infection [38]. Abaecin is a crucial immune effector involved in reactions against multiple parasites. Upregulation of several AMP genes, including abaecin and defensin-1, was also shown in other studies where honey bees were infested with mites [9,25]. Even though in our case Toll pathway was downregulated, the expression of abaecin can be also regulated by Imd pathway, through relish [39]. Second AMP, upregulated in whiteeyed pupae was defensin-1. Latest report showed that DWV could contribute to defensin 1 upregulation [25] therefore we suggest that defensin 1 could be modulated by DWV and different organisms that we may not have diagnosed, like bacteria, fungi, viruses and microsporidia. Defensin 1 is involved in the formation of social immunity, as peptide in royal jelly, by forming channels in plasma membrane that reduce integrity and permeability of cytoplasmic membrane in pathogenic organisms [38]. It appears that antimicrobial peptides, like Defensin 1, are able to protect honey bees even in early developmental stages and act as part of individual immunity [40]. Since defensin-1 is not regulated through Imd pathway [39] and in our case Toll pathway is not activated either, our results suggest that defensin-1 could be regulated through another pathway. In all the treatments down-regulation of iap-2 gene was noticed. Inhibitor of apoptosis 2 is known to be important for sustained response in the Imd pathway and its inactivation resulted in impaired microbial resistance in Drosophila [41]. Interestingly, the only thiamethoxam-dependent upregulation was detected for Dredd, a cytoplasmic caspase, involved in activation of Imd pathway, through cleavage of Relish [42].
In both groups of brown-eyed pupae, Varroa infested and simultaneously infested and thiamethoxam treated, the majority of genes became upregulated. In these two groups the cascades of genes involved in all four pathways, Toll, JAK/STAT, Imd and JNK, were upregulated ( Fig  1). It was reported [9,26,42] that peptidoglycan recognition proteins (PGRPs) are associated with mite parasitism and other pathogen infections. In our study, Varroa parasitism increased expression of PGRP-SC 4300, PGRP-LC 710 and domless. Activation of these three genes leads to activation of genes involved in all four pathways including effectors genes, defensin-2 and lysozyme-2 [16]. Contrary to earlier studies [19,20,23] our data show little evidence for reduced expression of immune-related genes in Varroa infested brown-eyed pupae. The only statistically significant downregulation was detected for two genes, spaetzle and PPOact. Prophenoloxidases (PPOact) are key enzymes responsible for activation of melanin synthesis, important in defense against microbes and wounding. Expression of this gene plays a crucial role in inhibition of parasite development in insects [9,43]. Spaetzle is known as a trigger for Toll pathway resulting in transcription of effectors (AMPs, PPOact and lysozyme) [26]. In brown-eyed pupae, thiamethoxam alone had only a minor effect on immune-related genes with no statistical significance.
The impact of in larval stage ingested thiamethoxam on Varroa infested white-eyed and brown-eyed pupae, had a significant effect that could be seen on heatmap (Fig 1) only when the mean expression level of each gene between Varroa infested and Varroa infested-thiamethoxam treated groups were compared (Fig 2). Normally, each treated group only to control group is compared. In infested white-eyed pupae, consumption of thiamethoxam leads to even stronger downregulation of spaetzle, abaecin and defensin-1, but it also caused higher expression of gene lysozyme-2, toll, and dorsal-1 when compared to group infested with Varroa (Figs 2 and 3B). It can be presumed that larval exposure to thiamethoxam deteriorates defense against Varroa mites and Varroa-mite associated viruses in pupal development stage. In brown-eyed pupae consumption of thiamethoxam to Varroa infested group leads to statistically significant increased expression of spaetzle, hopscotch, basket, and PPOact genes (Figs 2 and 3B). The only gene on which thiamethoxam had no effect in both white-and brown-eyed infested pupae was defensin-2. Expression of defensin-2 is induced by a pathogenic factor, and plays a role in individual immunity. Thiamethoxam caused more decline in gene expression on infested white-eyed pupae than in brown-eyed pupae that develop 3-4 days later.
During the infestation Varroa mites can vector viruses and other pathogens within and between colonies. Moreover 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 [22]. Comparable to other studies [9,16,21,44,45] a positive correlation between DWV and Varroa mites was found. We demonstrated that DWV RNA loads correlate with Varroa infestation in honey bee pupae (S1 Fig). Therefore, when talking about Varroa induced gene expression modifications, it would be important to investigate DWV impact on gene expression patterns [46][47][48].
Our data showed evidence for both down-and up-regulation of immune-related genes between two developmental stages of honey bees. Some of our results support studies reporting downregulation [18][19][20][21][22] other upregulation [9,16,25] of immune-related genes. It is difficult to explain such differential response in immune-related genes in two different pupal stages, however there is a link between the pupal development and immune response, a process of melanization. During the pupal period the intensity of melanization is rising with the development, in brown-eyed pupae this process being more intense than in white-eyed pupae. In immunity melanization is an immune effector mechanism involved in the killing different pathogens which involves a series of reactions to form a layer of melanin that surrounds and sequesters an invading pathogen [49]. Our results clearly show a difference in gene expression pattern between white-eyed and brown-eyed pupae, where significant change from downregulation to upregulation of almost all immune-related genes could be seen in the two to three days transformation. This study represents novel information about immune response to combined effect of thiamethoxam and Varroa parasitism. However, honey bee immune response depends also on post-transcriptional regulation as well as on post-translational regulation of gene expression and small differences in expression could reflect subtle modulation by a large number of factors acting in cascade [20]. Additionally, it is known that expression of genes can differ across developmental stages in different organisms [50]. Therefore, to paint a clear picture of pathogen and pesticide effects on honey bees each developmental stage should be analyzed.
In conclusion, our data provide novel insight into the genetic response of developing honey bees to exposure to thiamethoxam and Varroa parasitism. We demonstrated that Varroainduced activation of four known immune pathways does not occur earlier than in the browneyed pupal stage, and when thiamethoxam was present the significant alternation in expression of only 1/3 of analyzed genes was noticed. Thiamethoxam alone at concentrations within the  detected in hives had down-regulatory effects on immune-related genes in white-eyed pupae while in brown-eyed pupae its effect is minimal. In larval stage consumed thiamethoxam has an important and variable influence on immune related gene expression in Varroa-infested white-eyed and brown-eyed pupae. In white-eyed pupae down-regulation of several immune-related genes suggests increased thiamethoxam immune-toxicity after Varroa-infestation. In thiamethoxam intoxicated brown-eyed pupae up-regulation of major immune related genes indicate increase of immune activity after Varroa infestation. Our data indicates that thiamethoxam can contribute to Varroa harmful effects on honey bee pupae.