Brain metabolomic profiling of eastern honey bee (Apis cerana) infested with the mite Varroa destructor

The mite Varroa destructor is currently the greatest threat to apiculture as it is causing a global decrease in honey bee colonies. However, it rarely causes serious damage to its native hosts, the eastern honey bees Apis cerana. To better understand the mechanism of resistance of A. cerana against the V. destructor mite, we profiled the metabolic changes that occur in the honey bee brain during V. destructor infestation. Brain samples were collected from infested and control honey bees and then measured using an untargeted liquid chromatography-tandem mass spectrometry (LC-MS/MS)-based global metabolomics method, in which 7918 and 7462 ions in ESI+ and ESI- mode, respectively, were successfully identified. Multivariate statistical analyses were applied, and 64 dysregulated metabolites, including fatty acids, amino acids, carboxylic acid, and phospholipids, amongst others, were identified. Pathway analysis further revealed that linoleic acid metabolism; propanoate metabolism; and glycine, serine, and threonine metabolism were acutely perturbed. The data obtained in this study offer insight into the defense mechanisms of A. cerana against V. destructor mites and provide a better method for understanding the synergistic effects of parasitism on honey bee colonies.


Introduction
Honey bees provide pollination services to a diverse array of agricultural crop plants, which is a highly valued resource around the world [1]. The recent sharp decline in honey bee populations has caused a global crisis as the honey bee supply cannot keep up with the increase in agricultural demands [2]. Possible causes of the declines include exposure to certain pesticides, diseases, parasites, and even environmental deterioration [3]. Worth mentioning, the parasitic mite, Varroa destructor-originally confined to the eastern honey bee, Apis cerana-has become the most detrimental parasite of the western honey bee, Apis mellifera, and currently is considered as the major threat to apiculture [2]. PLOS  The factors by which V. destructor cause honey bee declines have been widely studied and have been attributed to two main causes: 1) female mites and their offspring, which feed on bee hemolymph for reproduction, weakening the honey bee's immune system, and 2) mites, which serve as an active vector of pathogenic viruses, and make the bees become more vulnerable to other secondary pathogens [4][5][6]. Previous studies have shown that all honey bee colonies are infested with Varroa mites. Unless steps are taken to reduce mite levels, colonies usually die within six months to two years. Meanwhile, infested bee species often exhibit dwindling populations and symptoms of viral and brood diseases [7].
Colonies of the Asian honey bee A. cerana are more resistant to V. destructor and suffer less damage than the western honey bee A. mellifera. Several effective defense strategies have been implicated such as grooming and hygienic behaviors [8]. Grooming behavior is regarded as the most important behavioral trait contributing to the defense against the parasitic mite. It consists of self-cleaning, a grooming dance, nestmate cleaning, and group cleaning [9]. Understanding the resistance mechanisms of the eastern honey bees against Varroa mites will help in breeding hybrids or other bees that are resistant to Varroa destructor. Dopamine, an important neurotransmitter in the brain, produces significant shifts in honey bee grooming behaviors and functions in a dose-and time-dependent manner [10]. Hygienic behavior is affected by olfactory cues such as by the odor of diseased, parasitized, or dead broods. The neuromodulator octopamine has been shown not only to enhance the response of bees to olfactory stimuli but also to play a role in biasing the nervous system for behavioral shaping [11]. Tyramine acts as the biosynthetic intermediate precursor for octopamine, which is capable of being released from neurons and influencing the honey bee's locomotor activity [12]. Taken together, these findings suggest that behavioral changes may be caused by neuropathological and neurophysiological effects on the honey bee's central nervous system.
Metabonomics is defined as the comprehensive quantitative and qualitative analysis of all metabolites in cells, tissues, or biofluids, and is as an emerging "-omics" method that provides surprisingly detailed insights into various metabolic processes under certain physiological or pathological conditions [13,14]. In the present study, a metabolomic analysis based on LC-Q-TOF MS/MS was first used to gain a better understanding of metabolic profile differences of A. cerana worker bee's brains before and after V. destructor infestation. The differentially expressed metabolites identified in this study will further increase our understanding of the resistance mechanisms of eastern bees against Varroa mites. Meanwhile, the data obtained in this study will also provide information likely to be critical for the development of rationally designed V. destructor prevention methods and for the breeding of new mite-tolerant honey bees.

Grooming behavior at the individual level
Mite resistance is considered a typical trait of eastern honey bees, and grooming is an important mechanism of protection against mite invasion. To determine grooming behavior changes during V. destructor infestation, both infested and non-infested bees were observed in a small modified Plexiglas 1 cage for 5 min, and individual grooming behaviors were recorded. As shown in Fig 1, the number of swiping behavior was extremely significantly higher in the V. destructor infection group compared to the control (p <0.001).

Analysis of brain metabolic patterns from infested and non-infested bees
To assess the capability of the LC-MS/MS-based metabolomic approach to differentiate V. destructor-infested bees from controls, we first analyzed all total ion chromatograms of honey bee brains. As shown in S1 Fig, the total ion chromatograms exhibited a stable retention time without obvious peaks' drifts. There were 8568 and 8210 ions identified in each sample profile under the ESI+ and ESI-mode, respectively. Nine QC samples were run for bee brain samples throughout the analysis. As shown in Fig 2, QC samples clustered together in the principal component analysis (PCA), which indicated that the LC-MS/MS system exhibited high stability and reproducibility. After removing low-quality ions (relative standard deviation (RSD) > 30%), we identified 7918 and 7462 ions in each sample analyzed in ESI+ or ESI-mode, respectively. To better illustrate the metabolic variations in honey bee brain after infestation by V. destructor, multivariate PCA and partial least squares discriminant analyses (PLS-DA) were applied to process the data. The plots of PCA scores showed no obvious separation between the infected group and the non-infected group in both ESI+ and ESI-modes (S2 Fig). Subsequently, the PLS-DA models used in further multivariate analyses exhibited good separation and are shown in Fig 3. Pairwise comparisons of abundances revealed ions in the infected group that were statistically differentially expressed from those in the control (corrected p-value < 0.05, FC > 1.2 or < 0.8). We constructed heat maps (commonly used for unsupervised clustering) based on these significantly different metabolites, which did not exhibit obvious segregation. In the ESI + mode, 130 ions differentiated the infection group from the controls (Fig 4A), and in the ESImode, 68 differentially expressed ions were detected ( Fig 4B).

Metabolic pathway analysis
We found putative identifications for 64 ion masses (all positive and negative mode ions) by searching against the mass-based HMDB database. Table 1 present the putative metabolites detected in ESI+ and ESI-mode, respectively. In the ESI+ mode, most identified metabolites were down-regulated. Among the down-regulated metabolites, iriomoteolide-1a was the least abundant metabolite (FC = 0.16). Meanwhile, more than half of the metabolites identified in the ESI-mode were less abundant. Metabolites like nizatidine, 2-propylglutaric acid, chlorphenesin, glycine, fonofos, and GMP were more abundant.
To identify the metabolic pathways disturbed during V. destructor infestation, the differential metabolites were subjected to KEGG pathway analysis and detailed related pathways were constructed using the reference map deposited in the KEGG database. Of the 25 differentially expressed metabolites detected in ESI-mode, 13 were mapped to the 12 KEGG metabolic pathways. Three metabolites were up-regulated, including GMP, glycine, and 2-propylglutaric acid. Further, of the 39 altered metabolites in ESI+ mode, 10 were mapped to the 10 KEGG metabolic pathways. The numbers of up-and down-regulated metabolites were six and four, respectively. Additionally, as shown in Fig 5A, an integrative view of the metabolic changes was prepared based on our findings. Furthermore, altered pathways based on all metabolites detected in both ion modes were analyzed with MetaboAnalyst 3.0 [15]. The impact value and-log (P) with metabolic pathway analysis (MetPA) were carried out to evaluate the importance of the pathways affected during V. destructor infestation. Among these altered pathways, we identified "linoleic acid metabolism", "propanoate metabolism", and "glycine, serine, and threonine metabolism" as three metabolic pathways of interest as all of them exhibited lower p-values and greater pathway impact ( Fig 5B).

Discussion
Tissue analysis is perhaps the most powerful approach for studying localized and specific responses to stimuli and pathogenesis and yields explicit biochemical information about the mechanisms of disease. LC-MS-based metabolomic technology has been widely used in the analysis of interactions between pathogens and their hosts [16][17][18]. The development of Varroa populations in honey bee brood cells often induces the death of a colony [19]. Parasitism by the mite causes Asian honey bees to exhibit a series of cleaning behaviors, especially grooming behavior, which can rapidly and effectively kill and remove the mite from the bee hive. In the score plot, each data point represents one bee brain sample, and the distance between points indicates the similarity between samples. x-and y-axes represent PC1 and PC2, respectively. https://doi.org/10.1371/journal.pone.0175573.g003 Brain metabolomic changes of honeybee infested with mite PLOS ONE | https://doi.org/10.1371/journal.pone.0175573 April 12, 2017 However, elucidating the molecular mechanisms underlying this hygienic behavior has been a persistent challenge for clinical and basic researchers. In the present study, LC-MS/MS based on a non-targeted metabolomic method was used to identify the characteristics of the metabolic profile of eastern worker honey bee brains infested with V. destructor. Our results provide important clues about the molecular metabolic network in the brain during mite infestation, especially the involvement of lipid, carbohydrate, and amino acid metabolism in the nervous system.  Grooming behavior is a crucial strategy for reducing ectoparasites in honey bees. While many studies have explored the mechanisms underlying the hygienic behavior, the molecular bias and functional pathways remain poorly understood. Navajas et al. [20] pointed out that the genes involved in neuronal development, neuronal sensitivity, and olfaction determine the susceptibility of honey bees to Varroa species. Later, using RNA-Seq, Mondet et al. [21] showed that antennae play an important role in the modulation of hygiene behavior, which functions through the alteration of the expression of antennal genes. Another study performed using proteomic methods identified proteins in the antennae of which the expression levels may influence the outcome of mite infestation [22]. Nevertheless, the underlying physiological and neuronal mechanisms remain unclear. Our data from the present study indicate that infestation-specific patterns were acquired in the brain metabolomics profile and that they can be used to distinguish honey bees with mite infestation from the controls. In addition, 64 putative metabolites were identified. Furthermore, by searching these metabolites against the KEGG database, many were found to be involved in metabolic pathways, including lipid metabolism, nucleotide metabolism, and amino acid metabolism amongst others. Integrated network analyses of the metabolites that were differentially expressed during mite infestation yielded highly related signaling networks, strongly suggesting that involvement of these signaling pathways could be essential for the resistance of mite infestation. Further, the metabolic network was mapped using MetaboAnalyst 3.0, and, based on our data, three metabolic pathways-"linoleic acid metabolism", "propanoate metabolism", and "glycine, serine, and threonine metabolism" were found to be the most relevant, especially linoleic acid metabolism is significant affected (p<0.05). Linoleate is a polyunsaturatedfatty acid, the abundance of which is found decreased in our study. Umezawa et al. [23] showed that male senescence-resistant SAMR1 mouse exhibited higher locomotor activity, more active circadian rhythm and higher number of response to stimuli after they were fed food high in linoleate. The majority of the differentially expressed metabolites detected through both ESI modes decreased in abundance in the infested bees compared to control bees. In the ESI+ mode, iriomoteolide-1a, fenapanil, phytolaccoside B, satratoxin H, and CDP-ethanolamine were most down-regulated with a fold change < 0.5. CDP-ethanolamine is involved in glycerophospholipid metabolism and the CDP-ethanolamine pathway, which is responsible for the de novo synthesis of phosphatidylethanolamine [24]. In humans, phosphatidylethanolamines are found particularly in nervous tissue such as the white matter of brain, nerves, neural tissue, and in spinal cord, where they make up 45% of all phospholipids [25]. Meanwhile in the ESI-mode, another four metabolites exhibited two times the down-regulated expression levels, including lolitrem K, cucurbitacin A, (1'S)-averantin, and L-histidinal. L-histidinal, a biosynthetic precursor of histidine, is involved in histidine metabolism and the biosynthesis of secondary metabolites [26]. In addition to the most altered metabolites, there are a few others worthy of mentioning. For example, GMP is a precursor for cyclic GMP (cGMP) through the action of phosphodiesterase (PDE), and our results showed that the level of GMP increased with mite infestation, which may function through the involvement of purine metabolism and hormonal signaling [27,28]. A previous study showed that GMP along with other guanine-based purines exhibited important neruromodulatory function [29]. In rat, acute administration of GMP was able to decrease the levels of anxietyin classical behavioral tasks [30]. Sphingolipids are highly bioactive compounds involved in diverse cell processes including cell-cell interactions and cell proliferation, differentiation, and apoptosis [31]. In this study, sphingosyl phosphochine was up-regulated in the sphingolipid metabolism pathway, suggesting sphingolipids may be an intracellular second messenger to function. Moreover, the alternation of sphingosyl phosphochine associated with Varroa-defense behavior mechanisms has not been reported; thus, further studies are needed to elucidate these.
Rhythmic behaviors-such as walking, hygienic behavior, grooming, and flying-are controlled by central pattern generators [32]. Grooming behavior significantly increased in bees infested by mites, indicating that some substances may act on neurons involved in pattern generation. Mosapride citrate hydrate has been recognized as a 5-HT 4 agonist, which exhibited decreased expression levels in this study. 5-HT 4 is a glycosylated transmembrane protein that functions in both the peripheral and central nervous systems to modulate the release of various neurotransmitters [33]. Isopentenyl phosphate (IPP), a down-regulated metabolite in the mevalonate pathway, functions as a dual inhibitor for transient receptor potential (TRP) A1 and TRPV3. Bang et al. [34] showed that peripheral IPP could administer attenuated TRPA1 and TRPV3 agonist-specific acute pain behaviors. Securinine, exhibiting a decreased level in infested bees compared to control bees in our study, serves as a selective antagonist of GABA recognition sites on mammalian central neurons, and it can inhibit GABA binding to rat brain membranes [35].
In conclusion, changes in the metabolic profiles of eastern honey bee brains exposed to V. destructor infestation were observed by a high-throughput non-targeted metabolomics method. Sixty-four putative significantly changed metabolites were identified, and the major metabolite network was predicted by pathway analysis. The putative metabolites were involved in a variety of pathways related to amino acid metabolism, nucleotide metabolism, and lipid metabolism. Worthy to mention, previous studies showed that the tactile and olfactory aspects of the mite were closely related to the bee defense behavior. However, the current study doesn't distinguish between tactile and olfactory cues that may be important for eliciting grooming behavior in bees. This will be our next study focus. The differential metabolites and related pathways identified not only further our understanding of the metabolic responses of honey bees to Varroa parasitism but also provide important molecular clues involved in the signaling cascades regulating hygienic behaviors. Also worth mentioning is that A. mellifera is susceptible to the V. destructor mite, frequently resulting in colony disintegration. Our present work lays a solid foundation for broadened the understanding of western honey bee metabolic profiles after mite infestation.

Honey bee colonies
The A. cerana populations were reared at the Apicultural Research Institute, CAAS, Beijing, China. Three honey bee colonies with mated queens of the same identical age and similar colony strength were selected as the experimental colonies.

Mite collection
Adult V. destructor mites were collected from a highly infested A. mellifera colony. Firstly, 50 adult A. mellifera were collected from hive frames and then transferred into a wire cage (8.7×6.0 cm), which was covered with a 4.0 mm wire mesh screen. Next, the wire cage was placed in a sealed container where the bees and mites were anesthetized with CO 2 for 2 min, allowing the mites to fall off from the worker bees. The mites were collected from the bottom of the container, placed in a petri dish, and kept at room temperature in the laboratory. This mite collection method was referred to in Arechavaleta-Velasco et al.'s study [36].

Observation hive
The observation hive used in this study was designed according to the method described by Arechavaleta-Velasco et al. [36] with slight modifications. Briefly, the observation hive was built out of wood (45×16×3 cm), and one end was covered with a piece of wire mesh (16×10 cm) for ventilation. The rest of the hive was covered with a transparent clear plastic sheet 1 mm thick. In the center of the plastic sheet, a circular hole 9 cm in diameter was made, which allowed the Plexiglas 1 cage (observe cage) to be inserted into the hive.

Behavioral assays
Approximately 1000 worker bees (A. cerana) in a comb along with one queen from a colony were firstly placed in the observation hive, which was then sealed and taken to the laboratory [36]. The schematic of the observation hive is shown in S3 Fig. Grooming behavior is defined as swiping motions in the direction of the mite with the front two pairs of legs. Grooming behavior at the individual level was evaluated as previously described [37]. Ten to 11-day-old worker honey bees were used in this study, and bees at this age show strongest grooming response. A total of 300 bees were divided into infestation group and control group, each group consists of 10 replicates. Each replicate having 15 bees, and these bees were equally selected from the three colonies. Bees in infestation group were selected and anesthetized with CO 2 to allow for manipulation of the bees. Then each bee was introduced into a small (10 × 10 × 7 cm) modified Plexiglas 1 cage. The Plexiglas cage has small pores (diameter 0.3 cm) on its wall, through which the tested bees can not walk out. When a bee appeared to have fully recovered, a V. destructor mite was placed on her thorax using a fine brush. The observation assays were performed in a dark environment with a red light on. Grooming behaviors were recorded with a video camera in 5 minutes. As for the control group, use the same method except for infested Varroa destructor. A Student's t-test was conducted to test for statistical differences in grooming behavior between control and infested mites; p-values < 0.05 were deemed statistically significant.

Sample collection
Another 300 honey bees were collected from the same test colonies and divided into two groups: V. destructor-infested and control groups, 10 replicates in each group. Each replicate having 15 bees, and these bees were equally selected from the three colonies. Using the same method, observation was performed on bees grooming behaviors for 5 min before being flashfreezing in liquid nitrogen.

Brain dissection and metabolite extraction
Using a dissecting microscope, bee brains were dissected on dry ice to remove all traces of the optic lobe and then stored at −80˚C until further processing. A total of 25 mg of brain samples was homogenized in 80 μL of a precooled methanol/water (1:1) solvent and ground for 5 min. The mixture was centrifuged at 20 000 ×g for 10 min at 4˚C. Supernatant (200 μL) was transferred to a new 1.5 mL polypropylene tube and processed through vacuum freeze drying before liquid chromatography separation.

Data processing
For quantitative metabolomics, raw data files were uploaded into Progenesis QI software (version 2.1), within which data alignment, normalization, and peak picking were performed. Next, the data matrix was mean-centered, pareto-scaled, and analyzed by both principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA). The quality of the models was evaluated with the relevant R 2 and Q 2 as described elsewhere [38]. Meanwhile, permutation tests (1000 cycles) were conducted to assess the robustness of the PLS-DA model. Differential metabolites were selected when the statistically significant threshold of variable influence on projection (VIP) values obtained from the PLS-DA model was larger than 1.0. Meanwhile, differences were tested with a Student's t-test, and a corrected p-value (q-value) < 0.05 was deemed statistically significant. Log 2 fold change (FC) was used to show how these selected differential metabolites varied between groups. A heat map was generated using the MultiExperiment Viewer (MeV) v.4.9 software (http://www.tm4.org/mev.html) based on the abundance of differentially expressed metabolite data (log 2 -scaled) [39]. To check and confirm the putative differentially expressed metabolites, the HMDB and KEGG databases were used [40,41]. Candidate metabolites were confirmed by MS/MS scans for the characteristic ions and fragmentation patterns of the compound. A pathway analysis of the differentially expressed metabolites was performed based on the HMDB and KEGG databases. MetaboAnalyst (http://www.metaboanalyst.ca) is a comprehensive web application for metabolomic data analysis and interpretation [15]. For the pathway enrichment analysis, a hypergeometric test was used for the over representation analysis, and relative-betweenness centrality was used for the pathway topology analysis. In the MetaboAnalyst analysis, a p-value of 0.05 was set as the threshold for significance.