Resistance and Vulnerability of Honeybee (Apis mellifera) Gut Bacteria to Commonly Used Pesticides

Agricultural and apicultural practices expose honeybees to a range of pesticides that have the potential to negatively affect their physiology, neurobiology, and behavior. Accumulating evidence suggests that these effects extend to the honeybee gut microbiome, which serves important functions for honeybee health. Here we test the potential effects of the pesticides thiacloprid, acetamiprid, and oxalic acid on the gut microbiota of honeybees, first in direct in vitro inhibition assays and secondly in an in vivo caged bee experiment to test if exposure leads to gut microbiota community changes. We found that thiacloprid did not inhibit the honeybee core gut bacteria in vitro, nor did it affect overall community composition or richness in vivo. Acetamiprid did also not inhibit bacterial growth in vitro, but it did affect community structure within bees. The eight bacterial genera tested showed variable levels of susceptibility to oxalic acid in vitro. In vivo, treatment with this pesticide reduced amplicon sequence variant (ASV) richness and affected gut microbiome composition, with most marked impact on the common crop bacteria Lactobacillus kunkeei and the genus Bombella. We conducted network analyses which captured known associations between bacterial members and illustrated the sensitivity of the microbiome to environmental stressors. Our findings point to risks of honeybee exposure to oxalic acid, which has been deemed safe for use in treatment against Varroa mites in honeybee colonies, and we advocate for more extensive assessment of the long-term effects that it may have on honeybee health.


INTRODUCTION
The deleterious effects of pesticides on non-target organisms have been a critical area of study in the last decades (Mancini et al., 2019). The vital ecosystem services pollinators provide prioritize assessment of pesticide toxicity in this group, and evidence of pesticide harm to honeybees continues to accumulate (Decourtye et al., 2004a;Gill et al., 2012;Goulson et al., 2015). The western honeybee, Apis mellifera, maintains a wide distribution, generalist foraging behavior and pollination competences that make it one of the most important species of pollinators across the world (Hung et al., 2018). Specifically, honeybees play a keystone role in pollination of natural ecosystems and agricultural crops (Goulson et al., 2015;Hristov et al., 2020). However, honeybee populations have declined globally over the past decades, with a continued reduction in colony numbers recorded in the United States, northwestern Europe, and Russia between 1940 and 2010 (Potts et al., 2010; US Department of Agriculture (USDA), , 2021Brosi et al., 2017), and a striking average annual loss of 30% of colonies reported by northern hemisphere beekeepers over the past decade (Brosi et al., 2017;vanEngelsdorp et al., 2017). A combination of anthropogenic stressors appears to drive honeybee declines (Meeus et al., 2018), including reduced habitats (Naug, 2009;Breeze et al., 2014), reduced genetic diversity (Espregueira Themudo et al., 2020), use of antimicrobials in apiculture (Tian et al., 2012;Li et al., 2017;Raymann et al., 2017), and pesticide application (Boncristiani et al., 2012;Williamson et al., 2014;Johnson, 2015).
Agricultural and apicultural practices expose honeybees to a range of pesticides that can damage their physiology, neurobiology, and behavior, and ultimately may result in colony decline. Two pesticide groups with potentially detrimental effects to honeybee health are insecticides, which are employed within their natural environment, and acaricides, which are applied within colonies. Commonly used neonicotinoid insecticides inhibit nervous system function by antagonizing acetylcholine receptors (Casida and Durkin, 2013). In honeybees, sublethal exposure to nitroguanidine neonicotinoids can alter locomotion (Williamson et al., 2014;Charreton et al., 2015), memory and learning (Decourtye et al., 2004b;Dacher et al., 2005;Shi et al., 2019), consequentially impairing flight ability, navigation (Tosi et al., 2017), and the proboscis extension response to sucrose (El Hassani et al., 2008;Thany et al., 2015). Due to the adverse effects, nitroguanidine neonicotinoids have been banned in several countries (European commission, 2013;Dewar, 2017), while the less toxic cyanoguanidine neonicotinoids remain widely used (Iwasa et al., 2004;Shi et al., 2019). Although cyanoguanidine neonicotinoids, such as thiacloprid and acetamiprid, are considered safer, they can still affect honeybee behavior, memory, and immune functions, and can be lethal at high concentrations (Iwasa et al., 2004;Thany et al., 2015;Brandt et al., 2017;Liu et al., 2020). Honeybees are also frequently exposed to acidic acaricides that are used to treat parasitic Varroa mite infestations (Sammataro et al., 2008). Oxalic acid is the most widely used acidic acaricide (Nanetti et al., 2003;Brodschneider et al., 2019) and although its mode of action remains unknown, it is presumed to be safe for bees as it is naturally found in honey (Bogdanov et al., 2002;Moosbeckhofer et al., 2003). Honeybees tolerate treatment concentrations of 3.5% (Rademacher and Harz, 2006), but exposure to higher concentrations can increase mortality and induce behavioral changes, such as reduced nursing efforts or general inactivity (Schneider et al., 2012;Rademacher et al., 2017). Oxalic acid treatment may also alter bee physiology, reduce pH in the digestive tract and the hemolymph (Rademacher et al., 2017), and create permanent lesions in the digestive tract (Martín-Hernández et al., 2007), like necrotic cells (Gregorc and Smodiš Škerl, 2007).
Until recently, most research on the effect of pesticides on honeybee health focused on the direct effects on the bees, but accumulating evidence suggests that effects extend to the honeybee gut microbiome. The honeybee microbiome is simple and conserved, with 8-10 core bacterial taxa that are omnipresent, regardless of geographical origin Moran et al., 2012;Ellegaard and Engel, 2019; Figures 1A,B). The microbiome modulates immunity against pathogens, partakes in digestion of pollen and in the neutralization of toxins (Engel et al., 2012;Kwong and Moran, 2016;Kešnerová et al., 2017;Raymann and Moran, 2018), promotes host weight and health, and mediates hormonal signaling (Zheng et al., 2016). Adverse effects of certain agricultural compounds on honeybee gut microbes have been documented. The herbicide glyphosate can perturb the absolute and relative abundances of dominant bacterial community members and increase honeybee susceptibility to pathogens (Motta et al., 2018;Blot et al., 2019). Chronic exposure to the highly toxic nitroguanidine neonicotinoids (e.g., imidacloprid and thiamethoxam) in laboratory settings can induce changes in gut bacteria community composition of healthy honeybees . However, these results have not been found when bees are returned to the hive after exposure (Raymann et al., 2018b). Thiacloprid, a cyanoguanidine neonicotinoid, potentially invokes dysbiosis of the gut microbiome (Liu et al., 2020), and although lower field-level concentrations may not alter colony performance, they reduce immune expression against pathogens . The effects of other commonly used cyanoguanidine neonicotinoids, such as acetamiprid, have yet to be investigated. Pesticides and antimicrobials used in apiculture to mitigate parasite and pathogen infections can also impact the honeybee gut microbiota (Li et al., 2017;Raymann et al., 2017Raymann et al., , 2018a. Acaricides such as coumaphos and tau-fluvalinate can influence the structure of the bacterial community in the gut (Kakumanu et al., 2016). Oxalic acid has antibacterial activity, including against Lactobacillus strains isolated from honeybees (Diaz et al., 2019), yet inhibition of other key core taxa and alteration of the gut microbial community have not been evaluated. Thus, our knowledge remains sparse despite the potential adverse effects of commonly used pesticides on honeybee gut microbes.
In this study, we address the potential effects of the pesticides acetamiprid, thiacloprid, and oxalic acid on the gut microbiota associated with honeybees. These compounds are still used on a global scale (Adjlane et al., 2016;Amulen et al., 2017;US Environmental Protection Agency Office of Pesticide, 2020;Begna and Jung, 2021), although thiacloprid was recently withdrawn from approval in Europe (European Commission, 2020; European Food Safety Authority (EFSA), 2020). Oxalic acid is a widely used treatment against Varroa (Rademacher et al., 2017). These pesticides have been deemed safe for bees, but insight into how they may affect honeybee gut bacterial communities is necessary to obtain a more complete risk assessment. We first performed in vitro assays to elucidate the potential direct negative effects of these compounds on main core bacterial members of the honeybee gut microbiome ( Figure 1C). The honeybee gut is compartmentalized and includes the crop (purple), midgut (dark blue), and hindgut (soft blue), which is further divided into the ileum and rectum (created and reproduced with permission from Biorender.com). (B) The relative abundance of core microbial genera in the honeybee gut microbiota (based on our dataset). (C) Overview of the in vitro experiment. Lawns of eight core bacteria were used to test the effects of three pesticides inoculated on 1 cm filter paper discs (structures drawn in ChemDraw). Inhibition was scored as the presence/absence of a distinct zone of inhibition forming around the disc. (D) Schematic of the in vivo experiment. For each of the three colonies, 20 sub-colonies were established each with 25 bees exposed to either of the three pesticides dissolved in sugar water or sugar water alone (control). After 7 days of exposure, three bees per sub-colony were randomly picked, their guts dissected and extracted, and the community composition established using amplicon sequencing.
Subsequently, we performed an in vivo experiment to investigate if the indications from the in vitro experiment translated to community effects on the microbiome and changes in honeybee health ( Figure 1D).

In vitro Assay of Gut Bacteria Susceptibility to Pesticides
To test for effects of the acetamiprid, thiacloprid, and oxalic acid on the growth of honeybee gut microbes, we obtained 15 strains from eight core gut bacteria of the western honeybee (A. mellifera) ( Figure 1C; Moran et al., 2012;Bonilla-Rosso and Engel, 2018;Ellegaard and Engel, 2019). Specifically, we obtained two strains of Lactobacillus Firm-4 (codes ESL0291 and ESL0292), two strains of Lactobacillus Firm-5 (codes ESL0183 and ESL0184), two strains of Lactobacillus kunkeei (ESL0216 and ESL0219), two strains of Bifidobacterium asteroides (ESL0199 and ESL0200), two strains of Gilliamella apicola (ESL0171 and ESL0172), two strains of Snodgrassella alvi (ESL0145 and ESL0252), two strains of Bartonella apis (ESL0058 and ESL0059) and one strain of Frischella perrara (ESL0215). Lactobacillus and Bifidobacterium strains were maintained on de Man, Rogosa and Sharpe agar (MRSA) and cultured anaerobically; Gilliamella, Snodgrassella, Bartonella, and Frischella strains were maintained on Columbia Blood agar (CBA) +5% blood medium. The former three genera were kept at 5% CO 2 , while the latter was maintained under anaerobic conditions. All strains were kept at 34 • C (cf. Kešnerová et al., 2016Kešnerová et al., , 2017. For standardized growth, five biological replicates of each strain were grown in liquid media, de Man, Rogosa and Sharpe (MRS) for Lactobacillus and Bifidobacterium strains, Columbia Blood (CB) for Gilliamella, Snodgrassella, Bartonella, and Frischella strains), and once OD600 reached 0.5 (after 1 day), 100 µl of each culture was plated in solid media in two 90 mm petri dishes to generate a bacterial lawn. Then, adapting from Balouiri et al., 2016, once the inoculum dried on the plate, four 5 mm sterilized paper filter discs were placed on the surface of the plate. The paper disc received 5 µl of the compound solution, and 5 µl of sterile water in the case of the controls. The concentrations chosen for each of these compounds were based on their solubility in water: acetamiprid (4.2 mg/ml), thiacloprid (0.1 mg/ml), and oxalic acid (95 mg/ml) to ensure that we tested effects of maximum concentrations, although these are likely to be far higher than those encountered in the field. This trial revealed that the bacteria were only sensitive to oxalic acid, and we consequently tested its effects at concentrations 0.095, 0.95, 9.5 mg/ml, which are conceivably more common doses for microbes to be exposed to than the very high maximum dose Schneider et al., 2012). Each strain was grown either without (controls) or with (treatments) a compound present in replicates of five. Presence/absence of inhibition was scored 2 days after exposure ( Figure 1C).

In vivo Assessment of Pesticide Consumption on Gut Microbiota
Honeybee Collection and Sub-Colony Set-Up Honeybees (A. mellifera) from three different colonies were obtained from the Department of Agroecology, Aarhus University, Research Centre Flakkebjerg in Slagelse, Denmark on the 7th of September, 2020. The colonies were opened, and the queen was located. Once found, bees were sampled from frames without the queen so that the colony would keep producing brood. Drones were avoided. Adult workers were picked from the bottom of the frame in an effort to avoid foragers and older bees. This implies that we most likely predominantly sampled younger bees working within the hive, with more stable gut microbiomes than foragers (Jones et al., 2018). Although it would have been ideal to sample comparably aged bees across the three nests, the consistent responses of the three colonies (see section "Results") suggest that differences in age are unlikely to have impacted our findings. The bees were then placed in Styrofoam boxes for transport to the University of Copenhagen. The three queens from the colonies were half-sisters.
At the University of Copenhagen, bees were anesthetized with CO 2 so that they could be separated into experimental groups. The Styrofoam boxes were put in plastic bags that were filled with CO 2 and, once the bees were inactive, they were sorted into their respective experimental groups. Bees from each colony were divided in five replicates (25 bees per replicate) for each experimental group and housed in rectangular plastic boxes of 18 × 11 × 6 cm, similar to previous work (Evans et al., 2009;Huang et al., 2014; Figure 1D). Each of the boxes had sixteen 2 mm diameter holes in the lid for ventilation, and a larger hole was drilled in the lid to allow placing a 15 ml Falcon tube for feeding. The Falcon tube contained 12 ml 50% sucrose diet without pesticides for controls and with one of the three pesticides for treatment sub-colonies. Each Falcon tube had three holes in the bottom from which the bees were able to feed (Huang et al., 2014; Figure 1D).
Compound concentrations in the sugar water provided to treatment sub-colonies were based on previous work. Thany et al. (2015) tested 0.1 µg acetamiprid/bee and El Hassani et al. (2008) tested 0.5 µg/bee. We decided to follow the former, because it had the concentration that was closest to what we expect bees to encounter in the environment (134 ppb, cf. Mullin et al., 2010). For the calculations, the onetime dose was multiplied by the number of bees per sub-colony (25) and the number of days that the experiment lasted (7). For thiacloprid, Mullin et al. (2010) detected traces of thiacloprid in 2% of pollen samples, at levels of 7.8 PPB. Zaworra et al. (2019) registered an IC50 (inhibitory concentration for half the replicates) of 0.19 µg. Tison et al. (2017) found that concentrations of 0.5-50 µg/ml did not increase mortality, but the highest did alter memory, we therefore used the lowest concentration. Schneider et al. (2012) used a 3.5% oxalic acid treatment (3500 µg/ml), which they reported to be one of the most common concentration beekeepers use in beehives . Rademacher et al. (2017) used a concentration of 50 µg/bee, considering this to be more representative of what was found in bee colonies (acknowledging the probable accumulation within colonies). Therefore, we used the second value and calculated the concentration in the same manner as for acetamiprid. Based on this, we provided bees access to sugar water with a concentration of 3.5 µg/ml of acetamiprid, 0.05 µg/ml of thiacloprid, or 1750 µg/ml oxalic acid. Although concentrations used thus vary greatly between compounds, they were chosen to best reflect ecological relevance.
After sub-colonies were setup, they were kept at room temperature with dim light. The sub-colonies were checked daily to count the number of dead bees and to check if sugar water was consumed or not. We did not remove dead bees during the experiment to minimize disturbance and to reduce the risk of escapees. The approximate volume of sugar water consumed was recorded on days one through four, after which the volume decreased such that it was not easily discernable and thus could not be recorded. The Falcon tube was refilled as needed with each treatment. The experiment lasted for 7 days, after which dead and live bees were collected and frozen separately. The gut dissections following the experiment were performed only on the live bees.

Gut Dissection, DNA Extraction, and 16S rRNA Amplicon Sequencing
Gut dissections were performed in sterile conditions on up to three bees per sub-colony, for a total of 146 bees (n control = 36, n acetemiprid = 39, n thiacloprid = 39, and n oxalicacid = 32). The dissection area was sterilized with 70% ethanol and the dissection materials with 96% ethanol and a flame, and the dissections were performed in the presence of a Bunsen burner (cf. Carreck et al., 2013;. The bee cuticle was sterilized prior to dissection by soaking the bee in a 1% aqueous solution of bleach for 3 min and then rinsing it in sterile purified water for three times 30 s (Binetruy et al., 2019). The bee gut was obtained by gently pulling the sting out of the abdomen, and the whole gut came attached to it. The gut was placed in a screw top 2 ml Eppendorf tube with 50 µl sterile PBS and frozen at −20 • C for later DNA extraction. In addition to the experimental bees, we included 12 bees from each colony that had been frozen immediately after collection (termed field controls). These field controls turned out to differ in community composition from our experimental controls, so we excluded them from the main manuscript analyses, but include them in a set of Supplementary Text, Supplementary Figures 1-3, 5-7, and Supplementary Tables 1, 2, 9.
For DNA extractions, we used the Qiagen Blood and Tissue kit (Qiagen, Hilden, Germany), following the manufacturer's protocol, but modified based on . For the tissue lysis, we used sterile 7 mm stainless beads and added 250 µl of 0.1 mm sterile crystal beads, running the bead beater at 30 Hz twice for 1 and 2 min, respectively. A total of 180 µl ATL buffer and 20 µl proteinase K were added to each sample, vortexed and incubated at 56 • C for 3 h on a rotor. After incubation, 4 µl RNase A was added and the sample was incubated for 15-20 min. The remainder of the protocol was as per the manufacture's protocol. The elution volume was 100 µl of buffer AE, and it was passed through the column twice in a joint elute.
Diagnostic PCR was conducted to confirm the presence of sufficient bacterial DNA in 143 samples. We used primers for the V4 region of the 16S rRNA gene ( V4.SA504: 5 -AAT GATACGGCGACCACCGAGATCTACACCTGCGTGTTATGG TAATTGTGTGCCAGCMGCCGCGGTAA-3 and V4.SB711: 5 -CAAGCAGAAGACGGCATACGAGATTCAGCGTTAGTCA GTCAGCCGGACTACHVGGGTWTCTAAT-3 ). PCR conditions were 94 • C for 4 min, 35 cycles of 94 • C for 30 s, 56 • C for 30 s, and 72 • C for 30 s, and a final extension step of 74 • C for 4 min. Amplifications were confirmed on a 2% agarose gel. Three negative controls without guts and a positive control with a cellular mock community standard (Zymobiomics, Nordic BioSite ApS, Copenhagen) were included in the DNA extraction, and three additional negative control were added during library preparation and sequencing. DNA was sent to the University of Michigan's Microbiome Core 1 for paired-end 250 bp 16S rRNA Illumina MiSeq amplicon sequencing with V4.SA504 and V4.SB711 primers.

Statistical Analyses
All statistical analyses were performed in R studio v. 4.0.3 R Core Team (2020).
For the bacterial inhibition assay, a generalized linear mixedeffects model (GLMM) with family binomial was used to address the effect of oxalic acid concentration on bacterial inhibition. Concentration was included as a fixed effect and bacterial strains as a random slope. Model reduction with subsequent ANOVA comparison was performed to address the effect of concentration on inhibition.
To analyze the mortality data (Supplementary Table 3), we used Cox proportional hazard regression models and likelihood ratio (LR) test statistics, employing the R package survival and the function coxph (version 3.2-7; Therneau, 2021). Treatment was included as a fixed effect, where treatments were compared to the control, and colony was included as a stratified fixed effect. Additionally, we created models to assess colony-level responses to treatment with treatment a fixed effect. For all models, the proportional hazard assumptions and the Cox-Snell residuals were tested according to Mills (2012); our models met both assumptions.
We used a GLMM with family binomial to assess if the type of pesticide affected whether or not the bees consumed sugar water. Treatment, colony, and day were inserted as fixed effects and sub-colony as a random slope. We then assessed the average volume of pesticide or control solution consumed by each bee over days one through four by dividing the daily consumption of the sub-colony by the number of bees alive in the sub-colony that day. We then used a general linear model (LM) to test the fixed effects of treatment, day, and colony on the volume of control or pesticide solution consumed. The minimal model was determined with backwards model reduction and met assumptions of normality and homoscedasticity. In both models, subsequent model reduction and ANOVA comparisons were performed to address the significance of the fixed effects on consumption and volume consumed.
For the 16S rRNA gene amplicon sequencing analysis, we used the dada2 pipeline (v.1.12.1; Callahan et al., 2016) and performed downstream analysis with default parameters with the following adjustments: we set the truncLen parameter in filterAndTrim to c(240,220), trimleft (6) and maxEE to c(2,3). We obtained 610 amplicon sequence variants (ASVs) after quality filtering and removing chimeric reads. The merged sequences were assigned to taxonomic ranks using the assignTaxonomy function in the dada2 package in R using the SILVA database release 138.1 (Quast et al., 2013). Since negative controls included core honeybee gut microbiome taxa, an abundance threshold was chosen to include genera in our study, while avoiding contaminants (cf. Kešnerová et al., 2017;Raymann et al., 2017). Only genera with an abundance >0.08% across the whole dataset were retained for further analysis, resulting in 15 genera and 203 ASVs. By removing non-abundant genera, we reduce the presence of putative contaminants, while keeping 98.7% of the total number of original reads, including genera that are part of the honeybee gut core microbiota and present in low levels in negative controls. The cellular mock community validated the detection of all eight expected taxa and allowed for an estimation of the random error in our community abundances, averaging 0.07% of the abundance of any given taxon.
Richness was estimated as the number of ASVs using the function estimate_richness in the R package phyloseq (v.1.34.0; McMurdie and Holmes, 2013). To test for statistical differences in richness between treatments and colony origin, a LM was used, where the log of the observed number of ASVs per treatment was the response variable. Both treatment and colony were included as fixed effects, and the model included their interaction. Pairwise comparisons between treatments were preformed using Tukey honestly significant difference (HSD) test, after confirming that the model assumptions were met using shapiro.test (v.0.9-38; Royston, 1982) and bptest (Breusch and Pagan, 1979) from the lmtest package (v.0.9-38; Zeileis and Hothorn, 2002). Beta diversity metrics were calculated based on Bray-Curtis and Jaccard distances, using the vegdist function from the vegan package in R (v.2.5-7; Oksanen et al., 2018). Additionally, Unifrac distances were also calculated using the UniFrac function. The differences between controls and each treatment were separately tested using multivariate PERMANOVAs (adonis function in the vegan package), together with a colony fixed effect and the interaction between colony and treatment. ALDEx2 (v.1.22.0;Fernandes et al., 2013) was used to analyze differentially abundant genera between controls and each pesticide treatment group. The same was done to detect differentially abundant ASVs. In each case, we controlled for colony by adding it as a fixed effect.
Lastly, as previous work has reported syntrophic relationships among honeybee gut symbionts (e.g., Kwong and Moran, 2016;Kešnerová et al., 2017), we performed network analyses using Bray-Curtis distance in the plot_net function in phyloseq for each treatment separately to test if associations in the control group were maintained with different pesticide treatments. We performed the analysis across all colonies to secure a sample size that allowed for a robust analysis. Distances below 0.5 are reported as associations.

Only Oxalic Acid Inhibited Bacterial Growth in vitro
We evaluated the proportion of plates with honeybee gut bacteria that were inhibited across the specified concentrations of three pesticides and found that all strains were resistant to acetamiprid and thiacloprid, but sensitive to oxalic acid (Figure 2). B. apis and Lactobacillus Firm-5 were only inhibited at the highest concentration (95 mg/ml), which is unlikely to be encountered in the field. Bacteria susceptible to concentrations relevant to acaricide application and bioaccumulation included S. alvi, F. perrara and B. asteroides, which were inhibited at 9.5 mg/ml, and Lactobacillus Firm-4, L. kunkeei, and G. apicola, which were most sensitive and inhibited already at 0.95 mg/ml; although, the latter for only one of the five test plates (Figure 2). The concentration of oxalic acid had a significant effect on inhibition (overall GLMM: F 2 , 220 = 25.13; concentration: χ 2 = 151.8, df = 3, p < 0.0001).

Treatment Affects Microbiome Composition, Particularly After Oxalic Acid Exposure
Across the full dataset, the most abundant genera were Snodgrassella, Lactobacillus, Gilliamella, Commensalibacter, Frischella, Bombella, Bifidobacterium, and Bartonella (Figure 5 and Supplementary Table 9), accounting for 87.2% of all clean sequence reads, and consistent with previous findings (Engel et al., 2012;Bonilla-Rosso and Engel, 2018;Ellegaard and Engel, 2019). There was, however, extensive variation both between colonies and, to a lesser extent, between bees from the same colony (Figure 5). Hafnia-Obesumbacterium and Klebsiella, two environmental bacteria that are opportunistically associated with honeybees, were abundant in some samples, with an overall average abundance of 3.9 and 9.7%, respectively. However, the abundance of these opportunistic bacteria was not associated with pesticide treatment (Figure 5 and Supplementary Table 9). Their elevated levels, including in controls, may have been due to the presence of dead bees within the boxes, which could have acted as reservoirs for the transfer of opportunistic and pathogenic bacteria to live bees. Future work should thus consider removal of dead bees to prevent effects on microbiomes of focal individuals.
To better understand the community changes underlying differences in alpha and beta diversity, we assessed the differential abundance of genera ( Figure 6A) and ASVs ( Figure 6B) between respective treatments and controls. We did not find any differentially abundant genera nor ASVs between controls and acetamiprid or thiacloprid treatments. However, six genera and seven ASVs differed in relative abundance between controls and oxalic acid treated bees. At the genus level, the relative abundance of Bombella was significantly reduced after oxalicacid treatment (ALDEx2 test: t = −6.054, p < 0.0001), while Gilliamella (t = 8.989, p < 0.0001), Frischella (t = 4.030, p = 0.0023), Lactobacillus (t = 7.986, p < 0.0001), Bifidobacterium (t = 6.912, p < 0.0001), and Snodgrassella (t = 3.942, p = 0.0033) all increased in relative abundance in response to the oxalic acid treatment (Figure 6A). At the ASV level, Bombella intestini (t = −8.326, p < 0.0001) and L. kunkeei (t = −10.56, p < 0.0001; Supplementary Figure 7) were negatively affected by oxalic acid. Conversely, a Lactobacillus Firm-5 ASV (t = 7.505, p < 0.0001), the most abundant Gilliamella ASV (t = 7.379, p < 0.0001), and three Bifidobacterium ASVs (t = 4.543, p = 0.0157; t = 4.210, p = 0.0251; and t = 4.067, p = 0.0348) were all relatively more abundant in oxalic acid treated bees than in controls ( Figure 6B). These changes in relative abundance of genera and ASVs reflect altered community composition. However, differences in the absolute abundances would require quantification of the 16S rRNA gene to estimate bacterial load, and should be considered in future studies.

Pesticide Treatment Affects Network Relationships Between Honeybee Gut Microbes
To elucidate potential syntrophic relationships among the identified genera in our dataset, and to assess whether these were affected by pesticide treatment, we performed network analyses for controls and each treatment, pooling data from all colonies (Figure 7). Overall, we observed flexibility in ecological networks of pesticide-treated microbiomes, but some positive associations were shared with controls (Figure 7). All positive interactions in controls were present in the acetamiprid treatment (Figure 7), but here we also saw associations between Commensalibacter and Lactobacillus, and between Snodgrassella and Frischella and Bartonella. In oxalic acid-treated honeybees (Figure 7), the association between Bombella and Snodgrassella present in controls disappeared, likely due to Bombella sensitivity to this pesticide. Additionally, Snodgrassella was linked to Lactobacillus, which potentially fills Bombella's niche of aerobic respiration. Under thiacloprid treatment (Figure 7), it is noteworthy that we could not replicate the Bifidobacterium-Frischella association, while a Lactobacillus-Snodgrassella association emerged. Finally, associations between Lactobacillus and other core members were surprisingly lacking in the controls, despite ample evidence of syntrophic relationships (Kwong and Moran, 2016). The lack of association may be due to the variable niches and interbacterial interactions of L. Firm 4, L. Firm 5, and L. kunkeei (Bonilla-Rosso and Engel, 2018) that lead to variable levels of Lactobacillus and inconsistent positive correlation with other bacteria. Nevertheless, overall changes in interbacterial networks point toward the presence of either a flexible or a fragile ecological network in the honeybee gut microbiome.

DISCUSSION
We examined the effect that three commonly used pesticides have on microbial symbionts of honeybees and found evidence for both in vitro and in vivo effects of oxalic acid exposure, but no effects of acetamiprid or thiacloprid exposure in the in vitro assay and far less marked effects in vivo. It is apparent that in vitro testing alone is insufficient to deduce effects on honeybee gut microbiomes and ultimately honeybee health. Although we were limited from quantifying the amount of liquid that bees consumed within sub-colonies late in the experiment, when fewer bees were present, all sub-colonies appeared to FIGURE 6 | Oxalic acid treatment changes microbial relative abundance in vivo. (A) Relative abundance of core microbes at the genus level between controls and different pesticide treatments. (B) Relative abundance of differentially abundant ASVs between control and pesticide treatments. Asterisks indicate significant changes between control and oxalic acid treatments assessed in two independent ALDEx2 bivariate models with pesticide treatment and colony as fixed effects (*p < 0.05, **p < 0.01, ***p < 0.001). Acetamiprid and thiacloprid are not significantly different from controls. Boxplots represent the first and third quartiles of the relative abundance, the horizontal line represents the median, whiskers extend 1.5 interquartile ranges and dots represent outliers.
consume sugar water, indicating exposure. Oxalic acid was the only pesticide for which we observed reduced number of days with consumption, consistent with previous suggestions that palatability may cause oxalic acid avoidance (Nanetti et al., 2015). Its consumption is nevertheless corroborated by the observed effects on gut microbial communities, with reduced microbial richness, reduced abundance of key microbiome taxa, and altered interbacterial relationships. We did not see consistently higher mortality induced by the presence of either of the three compounds compared to controls, so longer-term exposures may be needed to elucidate any such effects.
It should be noted that we observed very high mortality across all experimental groups, suggesting that caution should be taken when interpreting beyond laboratory experimental group comparisons. This high mortality was likely due to a combination of factors. The low alpha diversity in field control samples (see Supplementary Text) may suggest that the colonies were not in optimal health condition at the time of sampling (September). However, the alpha diversity and microbiome composition is known to change as the bees transition to the winter season, but the impact of this change on bee health remains unclear (Ludvigsen et al., 2015;Bleau et al., 2020;Kešnerová et al., 2020). Furthermore, although we conceivably avoided foragers (older bees) by sampling from the bottom of the hive frames, we did not strictly age control and may thus have included older bees with a higher risk of dying during the experiment. However, two other factors have conceivably played a more important role in governing mortality. First, honeybees are sensitive to stress (Klein et al., 2017), and transport and the use of CO 2 to anesthetize bees may have negative impacts. Secondly, to minimize disturbance and the risk of escapees we did not remove dead bees during the experiment, and it is likely that this has increased risks of cross infections with opportunistic pathogens from dead to live bees (see results on the elevated levels of opportunistic pathogens in control bees). Although such high mortality is not optimal, the consistent patterns observed in the microbiome analyses, and the impacts we see on specific core gut microbes, warrant meaningful comparisons and reflect biologically relevant effects.  (Anderson et al., 2013;Kwong and Moran, 2016;Butler et al., 2013;Kešnerová et al., 2017;Bonilla-Rosso and Engel, 2018).

Do in vitro Observations Predict in vivo Effects?
We did not observe in vitro inhibition by thiacloprid or acetamiprid of any of the tested bacteria, and the communities treated with these pesticides in vivo were only slightly altered compared to controls. Oxalic acid inhibition in vitro was consistent with greater alteration of the gut microbiomes in vivo, but the bacteria affected in vivo differed from those affected in vitro. This, and its comparably wider use as an approved pesticide, lead us to focus our discussion on potential impacts of oxalic acid exposure.
The outcome of pesticide treatment on gut bacterial abundances in vivo should be the combined effect of bacterial sensitivity, direct exposure and interbacterial dependencies. Bacteria that are sensitive in vitro to 0.95-9.5 mg/ml, conceivably comparable to potential exposures in the field, appear to increase in relative abundance in vivo after treatment with oxalic acid compared to controls. This is likely due to exposure differences in vitro and within bee guts, and illustrates the limitations that plating experiments have for evaluating impacts of anthropogenic compounds on host-associated microbes. Exposure is inevitable in vitro, but this is not necessarily so in vivo. For example, biofilm formation helps protect from environmental changes, including chemical or antibiotic exposure (Singh et al., 2017). Sensitivity to oxalic acid in the in vitro experiment would predict effects on e.g., G. apicola, F. perrara, B. asteroides, and S. alvi in vivo. However, within bees, G. apicola interacts in a biofilm with S. alvi and engages in cross-feeding interactions (Engel et al., 2012;Kešnerová et al., 2017). Such close syntrophic interactions within multispecies biofilms are likely to reduce oxalic acid exposure or buffer its effects.
Another potential explanation for the reduced impact on most core bacteria is their location in the intestinal tract. Apart from L. kunkeei and B. apis, the two microbes depleted by oxalic acid in vivo, all the tested bacteria reside in the hindgut (Anderson et al., 2013;Moran, 2015;Kwong and Moran, 2016;Powell et al., 2018). The hindgut receives predigested material from the midgut, and oxalic acid may thus at least partly have been digested before reaching the hindgut . A combination of decreasing oxalic acid concentrations along the bee gut and biofilm formation in the hindgut seems plausible to explain the effects on the honeybee gut microbes, suggesting that the most profound implications of exposure for colony health is conceivably effects on crop bacteria. In this context, we should note that increases in relative abundances are likely driven mainly by the loss of other bacteria, increasing the relative abundance of taxa in the microbiome that are unaffected by pesticide treatment. In order to determine if treatment affects the absolute abundance of bacteria, and to explore potential relationships between bacterial abundances and alpha diversity, bacterial load would need to be quantified.

Implications of Oxalic Acid Bacterial Inhibition on Colony Dynamics
In vitro inhibition of core honeybee microbes and in vivo depletion of crop members are potentially concerning for honeybee colony health. Bacteria move within colonies by trophallaxis between bees , which includes the exchange of crop contents, via bees picking up bacteria from hive material (Kwong and Moran, 2016), and through consumption of hindgut exudates (Powell et al., 2014). Oxalic acid treatment by trickling depends on the redistribution of the pesticide by the bees themselves (Schneider et al., 2012;Rademacher et al., 2017). While bacteria in most gut compartments may be somewhat protected from oxalic acid, we do observe a significant overall negative effect on microbiome richness, estimated as the number of ASVs bee guts contain. The regeneration of the microbiome after local extinctions in individual bees, or during inoculation of young bees, also depends on bacteria moving from hive to bee and between bees, implying that colony-level impacts from treatment with oxalic acid are likely.
Other studies have found oxalic acid inhibition in vitro of both honeybee microbes (Diaz et al., 2019) and plant pathogens (Kwak et al., 2016), and oxalic acid is present in honey and has been suggested to be responsible for its antimicrobial properties (Bogdanov et al., 2002;Nozal et al., 2003). While most strains appear resistant to oxalic acid in low concentrations, oxalic acid can persist in colonies for multiple weeks (Rademacher et al., 2017). If accumulation occurs after multiple oxalic acid treatments, the pesticide may have a stronger effect on honeybee microbes. Regular oxalic acid application to colonies may also select for resistance in opportunistic bacteria, Varroa, and honeybee mutualists and commensals. This may impair the function constitutive levels of oxalic acid have in honeybee colonies, ultimately requiring treatment with higher concentrations with impacts on colony health and production (Adjlane et al., 2016). As Varroa has a larger impact on honeybee health on the short-term than oxalic acid treatment, we advise exploring genetic resistance to Varroa in honeybees, as it may be a better long-term solution against infection than the application of oxalic acid (Conlon et al., 2019).

Implications of Pesticide Treatment on Honeybee Health
Oxalic acid impact on the honeybee gut microbiome could be direct or indirect. Indirect effects could be mediated by necrotic cell death in the midgut (Gregorc and Smodiš Škerl, 2007;Papežíková et al., 2017), lesions (Martín-Hernández et al., 2007), or pH changes (Rademacher et al., 2017). Since the affected microbes are present in the crop, rather than the midgut or hindgut, and the temporal scope of our experiment is restrictive, we can with some confidence rule out necrotic cell death and lesions for our experiment. This leaves pH changes as the most likely effector, which could have secondary effects on nutrient digestion and absorption.
Our network analyses replicate previously published data on interbacterial relationships (Kwong and Moran, 2016;Kešnerová et al., 2017;Ellegaard and Engel, 2019). Encouragingly, they revealed that interbacterial relationships are potentially flexible, at least on the short-term, with the caveat that our relationships are inferred by correlation, and hence may not reflect actual syntrophic or other types of relationships. Another interpretation is they are relatively fragile to environmental stressors. Crop bacteria are responsible for a portion of aerobic metabolism in the bee gut, therefore creating an anaerobic environment in lower portions of the gut. However, even after severe depletion of crop bacteria with oxalic acid, it seems that interbacterial relationships reassemble in such a way that Snodgrassella is now responsible for oxygen metabolism and depletion, leading to an association between Snodgrassella and fermentative Lactobacilli (Kešnerová et al., 2017). Snodgrassella's association with the facultative aerobe Gilliamella remains intact under this pesticide treatment, as do Gilliamella-Bifidobacterium and Bifidobacterium-Frischella associations, all of which exist in the lower portions of the honeybee gut. The main impacts on interbacterial relationships are therefore captured by the disappearance of Bombella and L. kunkeei, which affect other relationships, such as those involving Bartonella.
Lactobacillus kunkeei is one of the most abundant bacteria in the honeybee crop, and is also found in nectar, pollen, and honey (Anderson et al., 2013;Corby-Harris et al., 2014a;Bonilla-Rosso and Engel, 2018;Kwong and Moran, 2016). In vitro, we see that this species is relatively more resistant, potentially due to exposure of oxalic acid in honey (Bogdanov et al., 2002;Nozal et al., 2003). Lactobacillus presence increases brood and honey production (Alberoni et al., 2017), improves the honeybee immune response (Maruščáková et al., 2020), protects against foulbrood diseases (Vásquez et al., 2012) and some insect-associated strains can metabolize insecticides (De Almeida et al., 2017;Daisley et al., 2018). The strainlevel turnover between Lactobacillus strains we see in our experiment may have implications for colony health, as it is an important factor for honeybee microbiome function (Ellegaard et al., 2020). For example, L. kunkeei presence decreases larval mortality during Paenibacillus larvae infection and decreases Nosema infection prevalence in adults (Arredondo et al., 2018). Encouragingly, L. kunkeei can regrow surprisingly quickly after oxalic acid has been removed from the colony environment (Supplementary Figure 6). The other affected genus in our study, Bombella, has very specific niches within bees, including in nurse hypopharyngeal glands from which larvae are feed, the crop of nurse bees, and in royal jelly. Its location in the crop of nurses and in royal jelly implies potentially relevant roles in larvae and queen development (Corby-Harris et al., 2014b;Tarpy et al., 2015). Consistent with this assertion, it comprises a large fraction of the queen gut microbiome, possibly playing a role in queen nutrition and in modulating queen fertility, fecundity, and longevity (Anderson et al., 2018). This may involve protection from pathogens, as indicated by Bombella apis supplementation in colonies significantly reducing Nosema prevalence (Corby-Harris et al., 2016). Colony-level effects of oxalic acid treatment on Bombella must thus be assessed, as our results could have implications for honeybee immunity, longevity and fecundity.
Acetamiprid and thiacloprid treatment did not inhibit bacterial growth in vitro or impact the richness of honeybee gut microbiota compared to controls in vivo. Acetamiprid and thiacloprid treatments generate little change in the gut communities of treated bees, maintaining richness levels and core member relative diversity observed in the lab controls. Acetamiprid has not been investigated for its potential detrimental impact on honeybee gut microbiome, but chronic exposure of its fellow neonicotinoids, like thiamethoxam or imidacloprid, greatly alter intestinal communities of bees . Our in vivo findings contrast recent work by Liu et al. (2020) who described gut bacterial dysbiosis in honeybees exposed to thiacloprid in a dose dependent manner after one week of exposure; however, our experiment tested lower concentrations (0.05 mg/L), so this may explain the lower effect of thiacloprid. The bacterial communities recovered by day thirteen in Liu et al. (2020), likely due to anal trophallaxis between colony members compensating for the loss of honeybee gut members (Kwong and Moran, 2016) and could also be buffering the effect of thiacloprid in our experiment.

DATA AVAILABILITY STATEMENT
The 16S rRNA datasets generated in this study can be found in the SRA archive in GenBank under the BioProject PRJNA732842.