Honeybee Exposure to Veterinary Drugs: How Is the Gut Microbiota Affected?

ABSTRACT Several studies have outlined that a balanced gut microbiota offers metabolic and protective functions supporting honeybee health and performance. The present work contributes to increasing knowledge on the impact on the honeybee gut microbiota of the three most common veterinary drugs (oxytetracycline, sulfonamides, and tylosin). The study was designed with a semi-field approach in micro-hives containing about 500 honeybees. Micro-hives were located in an incubator during the day and moved outdoors in the late afternoon, considering the restrictions on the use of antibiotics in the open field but allowing a certain freedom to honeybees; 6 replicates were considered for each treatment. The absolute abundance of the major gut microbial taxa in newly eclosed individuals was studied with qPCR and next-generation sequencing. Antimicrobial resistance genes for the target antibiotics were also monitored using a qPCR approach. The results showed that the total amount of gut bacteria was not altered by antibiotic treatment, but qualitative variations were observed. Tylosin treatment determined a significant decrease of α- and β-diversity indices and a strong depletion of the rectum population (lactobacilli and bifidobacteria) while favoring the ileum microorganisms (Gilliamella, Snodgrassella, and Frischella spp.). Major changes were also observed in honeybees treated with sulfonamides, with a decrease in Bartonella and Frischella core taxa and an increase of Bombilactobacillus spp. and Snodgrassella spp. The present study also shows an important effect of tetracycline that is focused on specific taxa with minor impact on alfa and beta diversity. Monitoring of antibiotic resistance genes confirmed that honeybees represent a great reservoir of tetracycline resistance genes. Tetracycline and sulfonamides resistance genes tended to increase in the gut microbiota population upon antibiotic administration. IMPORTANCE This study investigates the impact of the three most widely used antibiotics in the beekeeping sector (oxytetracycline, tylosin, and sulfonamides) on the honeybee gut microbiota and on the spread of antibiotic resistance genes. The research represents an advance to the present literature, considering that the tylosin and sulfonamides effects on the gut microbiota have never been studied. Another original aspect lies in the experimental approach used, as the study looks at the impact of veterinary drugs and feed supplements 24 days after the beginning of the administration, in order to explore perturbations in newly eclosed honeybees, instead of the same treated honeybee generation. Moreover, the study was not performed with cage tests but in micro-hives, thus achieving conditions closer to real hives. The study reaches the conclusion that the most common veterinary drugs determine changes in some core microbiota members and that incidence of resistance genes for tetracycline and sulfonamides increases following antibiotic treatment.

micro-hives containing about 500 honeybees. Micro-hives were located in an incubator during the day and moved outdoors in the late afternoon, considering the restrictions on the use of antibiotics in the open field but allowing a certain freedom to honeybees; 6 replicates were considered for each treatment. The absolute abundance of the major gut microbial taxa in newly eclosed individuals was studied with qPCR and next-generation sequencing. Antimicrobial resistance genes for the target antibiotics were also monitored using a qPCR approach. The results showed that the total amount of gut bacteria was not altered by antibiotic treatment, but qualitative variations were observed. Tylosin treatment determined a significant decrease of aand b-diversity indices and a strong depletion of the rectum population (lactobacilli and bifidobacteria) while favoring the ileum microorganisms (Gilliamella, Snodgrassella, and Frischella spp.). Major changes were also observed in honeybees treated with sulfonamides, with a decrease in Bartonella and Frischella core taxa and an increase of Bombilactobacillus spp. and Snodgrassella spp. The present study also shows an important effect of tetracycline that is focused on specific taxa with minor impact on alfa and beta diversity. Monitoring of antibiotic resistance genes confirmed that honeybees represent a great reservoir of tetracycline resistance genes. Tetracycline and sulfonamides resistance genes tended to increase in the gut microbiota population upon antibiotic administration. IMPORTANCE This study investigates the impact of the three most widely used antibiotics in the beekeeping sector (oxytetracycline, tylosin, and sulfonamides) on the honeybee gut microbiota and on the spread of antibiotic resistance genes. The research represents an advance to the present literature, considering that the tylosin and sulfonamides effects on the gut microbiota have never been studied. Another original aspect lies in the experimental approach used, as the study looks at the impact of veterinary drugs and feed supplements 24 days after the beginning of the administration, in order to explore perturbations in newly eclosed honeybees, instead of the same treated honeybee generation. Moreover, the study was not performed with cage tests but in micro-hives, thus achieving conditions closer to real hives. The study reaches the conclusion that the most common veterinary drugs determine changes in some core microbiota members and that incidence of resistance genes for tetracycline and sulfonamides increases following antibiotic treatment.

RESULTS
General observations on the colony's status pre-and posttreatment. The trial involved bees treated with tetracycline (Pan-Terramicina, PT), sulfonamides (Sulfac, SUL) and tylosin (Tylan, TL), plus an untreated control (CTR); each experimental condition was tested with 6 replicates. Bees were sampled at T0 (experiment beginning) and T1 (24 days later). Moreover, the experiment relied on micro-hives managed with a semi-field approach due to national restriction on antibiotics.
Throughout the trial, no particular changes or deficiencies in the health status of the treated honeybees were observed. Only one micro-hive collapsed (PT_6) just after the experiment's end, presumably due to varroosis, whereas CTR_5, PT_1 and SUL_1 were found to be queenless at the end the experiment. Visual evaluation at the time of gut sampling highlighted a reddish coloration of the intestinal epithelium in the tylosin treatment group. Drought conditions in the second half of the experiment did not allow nectar harvest and consequently no weight increase was observed, despite the sugar syrup supplementation.
qPCR quantification of target microbial groups in the gut and resistance genes. The counts of Eubacteria (Fig. 1A) at the beginning and at the end of the experiment showed a significant decrease (0.65 log, P , 0.05) upon sulfonamide treatment (SUL_T0 versus SUL_T1). Other conditions did not show significant variation. Considering Bartonella spp. (Fig. 1B), only PT treatment highlighted a significant decrease between PT_T0 vs PT_T1 (0.76 log 16S rRNA copy number decrease, P , 0.05). Bifidobacterium spp. counts showed a general decrease in all experimental conditions. The reduction was significant in PT_T0 versus PT_T1 (0.58 log 16S rRNA copy number decrease, P , 0.01) and in TL_T0 versus TL_T1 (3.61 log 16S rRNA copy number decrease, P , 0.01) (Fig. 1C). Also, Bombilactobacillus spp. and Lactobacillus spp. showed a general decrease in all experimental conditions, which was significant only in the comparison of TL_T0 versus TL_T1 (P , 0.01), with a decrease of 2.89 and 1.71 log 16S rRNA copy numbers, respectively ( Fig. 1D and E).
Tetracycline resistance genes tetW and tetY increased significantly by 144% and 180%, respectively, (P , 0.01) comparing PT_T0 with PT_T1. Sulphonamides resistance Data are expressed as the log of 16S rRNA gene copies/intestine for Bartonella sp., Bifidobacterium spp., Bombilactobacillus sp., and Lactobacillus sp.; for Eubacteria, data are expressed as the log of 16S rRNA copies/intestine. Boxplots report minimum and maximum values, lower and upper quartile and median. CTR, no antibiotics control; PT, oxytetracycline; SUL, sulfonamides; TL, tylosin. genes sul1 and sul2 showed a significant increase (76.84% and 33.95%, respectively, comparing SUL_T1 with SUL_T0, P , 0.01), whereas sul3 could not be amplified at the different annealing temperatures tested (40 to 64°C). Tylosin resistance genes tlrB and tlrD did not show any significant variation in normalized data. The melting temperature (T m ) of the amplification products immediately after the last reaction cycle and the qPCR efficiency data are reported in Table 1.
Bee gut microbiota analysis via next-generation sequencing. A total of 48 samples (2 sampling times [T0 and T1], 4 experimental conditions [CTR, PT, SUL, and TL], 6 replicates for each condition, each replicate being a pool of 30 honeybee guts) were subjected to next-generation sequencing (NGS) analysis on an Illumina MiSeq platform. About 13.7 million raw reads were obtained from the sequencing. Of these, 9.1 million reads passed the quality control and the chimera check analysis, obtaining an average of 95,986 joint reads per sample. For statistical analysis, samples were rarefied at 48,400 reads, a value obtained with exclusion of one replicate (TL_T1_4) due to a particularly low coverage. The taxonomical assignment of the 47 samples produced 17,194 operational taxonomic units (OTUs) at 97% similarity based on the SILVA 132 database. The obtained NGS data on the whole data set are reported in Table 2, as well as absolute abundance at phyla, families, and genera levels per treatment and time. Figure 2A reports absolute abundance at genus level per replicate.
The a-diversity indices (Chao1, observed OTU, and PD whole tree) showed a significant decrease over time only in the tylosin-treated group (P , 0.01). The b-diversity indices, considering unweighted UniFrac, underlined statistically significant differences between the CTR and TL treatments. However, considering the abundance of taxa in the weighted UniFrac, not only TL treatment but also SUL treatment resulted in significant differences when compared to CTR.
The intestinal microbial taxa at the different taxonomic levels did not show any significant shift between the two sampling times (T0 and T1) in control bees. A summary of the significant changes, from phyla to species, for each antibiotic treatment over time is reported in Table 3.
Regarding tylosin treatment (comparing TL_T1 versus TL_T0), Proteobacteria doubled their abundance (P , 0.01). On the other hand, both Firmicutes and Actinobacteria significantly decreased (P , 0.01). Bifidobacteriaceae and Lactobacillaceae significantly decreased comparing TL_T1 and TL_T0 (P , 0.01), with percentage values that are consistent with those reported below at the genus level. Orbaceae significantly increased at T1 (168.63%, P , 0.01). Finally, the absolute abundance of Other_families significantly increased after TL treatment (1673%, P , 0.01). The Bifidobacterium spp. absolute abundance reduction after TL treatment was highly significant (P , 0.01), decreasing from 9.32% at T0 to 0.02% at T1 (Fig. 3B). In the same way, Bombilactobacillus spp. and Lactobacillus spp. decreased from 10.61% and 37.52% at T0 to 0.81% and 9.37% at T1 (P , 0.01) ( Fig. 3C and G), respectively. Moreover, the absolute abundance of Bartonella spp. and Gilliamella spp. strongly increased (P , 0.05) ( Fig. 3A and F). Other_genus species significantly increased from 1.84% to 12.40% (P , 0.01) (Fig. 3I). At species level, a significant decrease of six Lactobacillus species and also of unclassified Lactobacillus spp. was observed (P , 0.01), together with the decrease of B. mellis (P , 0.01), Bifidobacterium asteroides (P , 0.01), and Bifidobacterium indicum (P , 0.05). The Cramer V test showed a strong biological relevance in pairwise comparisons of TL_T1 versus TL _T0 and SUL_T1 versus SUL _T0 (Cramer  Taxon  T0_CTR  T1_CTR  T0_PT  T1_PT  T0_SUL  T1_SUL  T0_TL  Within the principal-component analysis (PCA) of the data set at species level, PC1 and PC2 together explained only 25% of the variability. However, the TL_T1 group is clearly separated from TL_T0 and also from the other treated samples at T1 (Fig. 4A), particularly along the PC1 axis. Orbaceae and thus Gilliamella spp. are associated with TL_T1, as also confirmed by statistical analysis (Fig. 4B and C). The graph also shows a clear separation of SUL_T0 and T1 along PC2.

DISCUSSION
In this work, we investigated the gut microbial community of honeybees after the administration of antibiotics (oxytetracycline, sulfonamides, and tylosin) against common bee diseases.
Total bacteria counts were not greatly affected by the antibiotic treatment, whereas the amount of some microbial groups varied significantly upon target antibiotic exposure.
Oxytetracycline is a broad-spectrum antibiotic currently used in the beekeeping sector (19,33). Recently, Raymann et al. (22,23) showed that the use of tetracycline strongly decreased the absolute abundance of 5 core gut genera in partially caged honeybees, in particular Bartonella, Bifidobacterium, Bombilactobacillus spp. (formerly known as Lactobacillus Firm-4), Lactobacillus, and Snodgrassella. In our study, the large increase of tetracycline-resistance genes in the gut bacteria upon antibiotic treatment is accompanied by the increase of some core members, two of which are significant (Gilliamella spp., in agreement with Raymann et al. [22], and Snodgrassella spp.). The abundance of other microbial genera, such as Bartonella and Bifidobacterium, decreased. Our results therefore show an important effect of tetracycline that is focused on specific taxa with minor impact on alfa and beta diversity. It is well known that honeybee gut commensal bacteria provide large reservoirs of tetracycline-resistance determinants (otr and tet genes) frequently acquired through large and/or long-term antibiotic exposure or from other habitats shared with animals and humans (34,35). Ludvigsen et al. (35) showed that honeybee gut symbionts, in particular Snodgrassella spp. and Gilliamella spp., can survive and proliferate thanks to tet determinants, and this further supports our results that show a significant increase of these two genera. Recently, Daisley et al. (36) found that the routine administration of oxytetracycline increases tetW and tetY abundance in the gut microbiota of adult workers and is associated with a depletion of the major symbiont taxa. The present study, therefore, confirms that honeybees may represent a reservoir of tetracycline resistance genes (Fig. 5). In addition, bees, with their daily activities (hive interaction, flying, flower visiting), have a preferred path to integrate their gut microbiota and mitigate the antibiotic damages upon tetracycline administration, as suggested by Daysley et al. (36) for Lactobacillus strains. Most of the published studies rely on caged or partially caged honeybees, which limits social behavior as well as interactions with environmental bacteria. In our study, an additional microbial source for the gut microbiota may derive by the reservoir of microbial inoculants within the hive structure (stored pollen, nectar, and wax), which may have contributed to the mitigation of tetracycline impact. Another mechanism that has to be considered in antibiotic resistance genes (ARGs) spreading regards the transmission through horizontal gene transfer (HGT) that has been documented in soil following application of manure containing antibiotic residues and also in human intestine. The HGT of antibiotic resistance genes may allow the increase of ARGs in the studied environment even without major perturbations of the microbiota (37-39). Sulphonamides (SUL) were widely used in the beekeeping sector from 1960 to 2000, but residues in honey are still found, thus showing they are still used in spite of the banning (40). Among the core genera found in the honeybee gut, Frischella and Bartonella spp. were significantly affected by SUL treatment, while Bombilactobacillus spp. and Snodgrassella spp. increased their counts. Frischella perrara has implications in immune priming in honeybees and in the induction of peptides with antimicrobial activity (41). The registered 3% reduction (with a final 1% abundance in T1) may have controversial implications. F. perrara reduction could be detrimental for the bee immune stimulation (41), on the other hand, this species has been reported as pathogenic because it causes scab formation in the pylorus (42), therefore Frischella reduction might also be positive. Bartonella spp. has been related to the recycling of nitrogenous waste products into amino acids and with the degradation of secondary plant metabolites (43). The reduction of more than 80% of this taxon could have implication in digestion functions and in the recovery of amino acids (43). However, it is evident that most of the core members are not affected by SUL treatment. This can again be a consequence of the increase of the sulfonamides-resistant population upon selection after sulfonamides exposure. Accordingly, Cenci-Goga et al. (44) found a high abundance of sulfonamides resistance genes (sul1 and sul2) in honeybees sampled in different Italian locations because of the high SUL spread in the environment.
Tylosin induced a remarkable change in some microbial taxa proportions, almost causing the depletion of the rectum population, in particular of lactobacilli and bifidobacteria, and favoring the hindgut population (mostly Gilliamella, but also Snodgrassella and Frischella). It is known that tylosin targets are mainly Gram-positive bacteria (45,46). Bifidobacterium, Bombilactobacillus, and Lactobacillus genera represented 99.99% of Bifidobacteriaceae and Lactobacillaceae family members that, overall, accounted for more than a half of the honeybee gut microbial community. They play an essential role in the transformation of various pollen coat-derived compounds, including flavonoids, phenolamides, and v-hydroxy acids (47), in addition to the digestion of complex sugars (48,49). Their rapid decrease may affect honeybee ability to metabolize specific compounds and consequently reduce nutrient availability. It is interesting to report that the macrolide antibiotic resistance genes tlrB and tlrD did not increase significantly in treated honeybees at T1, even if detected. This is probably due to the low occurrence of these ARGs in Bombilactobacillus, Lactobacillus, and Bifidobacterium honeybee strains, even if TL-resistant strains have been described in humans and swine (50,51). Tlr genes belong to the same resistance group as erm genes (erythromycin ribosome methylation), so that tlrB is also classified as erm32 and tlrD as ermN (52,53). The maintenance of tlr gene abundance may also be explained by their activity against other macrolide antibiotics with a broader spectrum of activity, including Gram-negative bacteria that survived the TL treatment. Indeed, Jackson et al. (54) found that erm genes can be activated after tylosin use. Several studies showed that environmental species, such as members of the Asaia, Apibacter, Apilactobacillus, Vagococcus, Pseudomonas, Parasaccharibacter, Citrobacter, Providencia, and Pantoea genera, often related with soil, pollen, and nectar (55,56), are present in the honeybee gut as minor groups (57)(58)(59). These non-core genera were found to increase at T1 upon treatments with SUL and TL. These microorganisms may promote the increase of the pool of ARGs due to their continuous exposure to antibiotics used in agriculture, such as the use of sewage from livestock as a soil amendment. Among these strains, Parasaccharibacter apium, recently reclassified as Bombella sp. by Smith et al. (60), is reported as a strong immune-stimulating strain in honeybees, also capable of counteracting Nosema sp. (61). Therefore, the non-core genera that are sporadically associated with honeybees might play a role in the immune stimulation or metabolic regulation of honeybees, despite their low abundance, and may increase upon antibiotic treatment.
Overall, the three assayed veterinary drugs seem not to influence the total amount of bacteria but rather the absolute abundance of several core and non-core taxa, causing a possible lack of metabolic functions related to the most susceptible bacterial species and strains. A long-term observation of the colony health status, also including the hive development and hive products (e.g., honey), will allow the understanding of the relationship between the altered microbial structure and the behavior and performance of honeybees.

MATERIALS AND METHODS
Experimental design. Due to the European and national laws restricting the use of antibiotics or other veterinary drugs in the open field, tests were conducted in semi-field conditions, i.e., in microhives incubated in a thermostatic chamber with a short flying time for honeybees. Honeybees employed in this study had not been treated with antibiotics for several generations (over 2 decades).
The micro-hives employed in the study were obtained as depicted in Fig. 6. A number of wax combs obtained from a fully populated and healthy bee colony were shaken on a box containing 72 new micro-combs (L 9.5 Â H 10.5 cm), causing the fall of thousands of honeybees on the provided microframes and populating them (this procedure is referred to as the "shook swarm" method). The queen was allowed to lay eggs for 3 days on approximately 1/3 of the total available micro-combs. Five days later, 24 experimental wooden micro hives (L 20 Â H 15 Â W 16 cm) were set up, each containing 3 micro-combs (a brood frame, a honey frame, and an empty comb). Each micro-hive contained approximately 500 honeybees with a mated queen. The obtained micro-hives constituted the experimental replicates (6 for each experimental condition). Moreover, every micro-hive was equipped with an anti-robbing entrance modification, forcing honeybees to cover an "S" path that discouraged the entrance of robber bees when the micro hives were placed outside.
Micro-hives were placed into an incubator with controlled temperature and humidity (29°C and relative humidity [RH] of 60), and equipped with a net allowing ventilation on the mini-hive bottom. The micro-hives were moved outside in the late afternoon (approximately from 5:30 p.m. to 8:30 p.m.) every second day in order to allow the bees to fly freely and defecate. The arrangement of the micro-hives outdoors in the experimental field always followed the same pattern to avoid disorientation and drift. Microhives were placed at minimum 2 m distance from each other, and in clusters of 3 units of the same experimental thesis, oriented in different directions, in an experimental forest well populated by trees.
At evening, micro-hives were closed and relocated to the lab incubator. Micro-hives were fed every 2 days with 30 ml 1:1 (wt/wt) sucrose solution, plus a dispenser containing 5 ml sterile water. The day of the antimicrobial treatment, honeybees were treated as described below. The developed experimental conditions were: TL, tylosin; PT, oxytetracycline; SUL, a mixture of sulfaquinoxaline and sulfadimethoxine; and CTR, the control with no antibiotic administration. Details on antibiotic use and concentrations are reported below.
The trial was carried out between July and August 2016, where two foraging options were available: blooming all through the trial, even if strongly limited by summer drought. The health status (adult honeybee population and brood size, honey reserves, core colony cohesion, symptoms of viral diseases, and varroa infestation) of honeybee micro-hives was periodically assessed, and variations annotated when relevant.
Treatment preparation, administration, and sampling. Antibiotics were administered according to available guidelines for each antibiotic (62)(63)(64). Details and concentrations of antibiotics are reported in Table 4. Bees were treated once a week for a total of three treatments with micro-hive feeders containing 30 ml of sugar syrup (1:1 wt/wt) mixed with the respective treatment. Finally, after the 3rd treatment (days 15 to 17), at least 50 emerging honeybees per replicate were marked on the thorax (65) with colored nail polish nontoxic to bees. Marked honeybees were sacrificed at day 24, at nurse stage (7 to 9 days post eclosure), and with a completely established gut microbiota (66). A pool of 30 bees per replicate (a total of 180 samples/experimental condition) was picked at the beginning of the experiment (T0) and after 24 days (T1).
DNA extraction and NGS sequencing. Obtained honeybee gut pools were well homogenized with a pestle, after which was added 1,400 ml of lysis solution containing 60 ml proteinase K per pool (20 mg/ ml concentration) and glass beads. Total destruction of gut epithelial tissues was obtained after 1 h incubation at 55°C. Only 1/4 of the resulting sludge (450 ml) was used for gut genomic DNA extraction with Quick-DNA Fecal and Soil Microbe kit (Zymo Research, California, USA). The 16S rRNA gene amplification and libraries prepared for Illumina MiSeq platform sequencing were carried out according to Alberoni et al. (67). Briefly, the V3-V4 was amplified with KAPA Hi-Fi PCR Master Mix (Roche, Monza, Italy) with a maximum of 25 cycles. PCR products were purified with AMPure magnetic beads (Beckman Coulter, Milan, Italy) and indexed with i7 and i5 Illumina adapters (Illumina, Milan, Italy). NGS sequencing was performed with the addition of 22% PhiX to the sample pool (Illumina, Milan, Italy). Bioinformatic analyses were performed with Qiime1, and representative operational taxonomic units (OTUs) were subjected to BLAST search against the most updated SILVA database release 132. OTUs with less than 0.1% abundance were discarded. The a-diversity was evaluated using Chao1, observed OTU, and PD whole tree metrics, whereas b-diversity was evaluated using both weighted and unweighted UniFrac.
Quantification of target microbial groups and resistance genes. The main microbial groups found in the honeybee rectum (Bartonella spp., Bifidobacterium spp., Bombilactobacillus spp., and Lactobacillus spp.), as well as total bacteria (Eubacteria) were quantified with qPCR (StepOne real-time PCR system, Applied Biosystems) according to Baffoni et al. (68,69). Briefly, standard curves were constructed using PCR products of the 16S rRNA gene for the target microbial genera. The PCR products were purified and serially diluted to obtain standards ranging from 10 4 to 10 8 gene copies. Amplification reactions were performed on a total volume of 20 ml using the Fast SYBR Green Master Mix (Applied Biosystems). Moreover, DNA quantity was standardized at 5 ng/ml of DNA. The specificity of the reaction given by the T m of the amplification products is reported in Table 1.
Absolute quantification of the target microbial groups was obtained by multiplying the absolute quantification data with the total extracted DNA and then divided by the gut number and the gene copy number (except for Eubacteria) (70,71). The data was expressed in log 16S rRNA copies/intestine (72). ARGs TetW, TetY, Sul1, Sul2, Sul3, TlrB, and TlrD (Fig. 5) were quantified according to Zhang et al. (73). The primers used are reported in Table 5. Raw data were corrected according to the total DNA quantification. The final absolute abundance of ARGs was normalized according to references 74 and 75 by dividing the total ARGs with the absolute abundance of total bacteria previously obtained.
Data adjustments and classification of microbial genera. Rarefied biom tables obtained from NGS bioinformatic analysis were used for further data adjustments, where the absolute abundance of each bacterial species was calculated according to Raymann et al. (22) by multiplying absolute abundance data to the corresponding qPCR total amount results, and then normalizing by the copy number of the 16S rRNA gene typical of each microbial genus. Moreover, species belonging to the Lactobacillus genus have been recently reclassified (76) but databases for NGS OTUs assignment were not yet updated with the new classification at the time of the bioinformatic analysis of the presented data. Therefore, the data set was manually adjusted according to Alberoni    Pro805-R plantarum to the new respective taxonomical classifications Apilactobacillus kunkeei and Lactoplantibacillus plantarum. Due to the concern that sequencing amplicon length (%470 bp) might not be enough to efficiently discriminate among species, manual curation was then used to validate by qPCR with Firm-4-and Firm-5-specific primers (30). The obtained data set was used for further graphical and statistical analyses on target genera and species. Compliance with ethical standards. This article does not contain any studies with human participants by any of the authors and experiments on animals were performed according to the Italian laws allowing experiments on arthropods without the need of an official ethical commission approval, unless cephalopods are used.
Statistical analysis. Statistical analysis for qPCR and NGS data (a-diversity and taxon analysis) was performed with the R software (78) according to Baffoni et al. (68). Analysis on data normality and homoscedasticity was performed and normal and homoscedastic data were analyzed with ANOVA; nonnormal homoscedastic data (with normal distribution of residuals) were analyzed with glm function, while data with high deviation from normality were analyzed with the nonparametric Kruskal-Wallis test coupled with the Dunn test. For b-diversity index, data resulting from QIIME statistical elaboration were reported. The software calculates the UniFrac distance (weighted and unweighted UniFrac) between all the pairs of samples in the data set to create a distance matrix. The statistical significance between groups was subsequently estimated using the Monte Carlo method with the Bonferroni correction.
Post hoc tests among different groups were carried out and Bonferroni's correction was applied. The post hoc test considered pairwise comparisons within each experimental condition, taking into consideration the impact of each treatment over time. Therefore, four comparisons for the semi-field trial and three comparisons for the in-field trial were considered. The control was considered as a further treatment to monitor and evaluate the normal gut microbial community evolution resulting from the interaction of honeybees with the environment. Graphs were generated with ggplot2, ggpubr, and Microsoft Excel. The biological relevance of experimental conditions, pairwise compared at their respective sampling time (T1 versus T0), was computed with Cramér's V (79) relying on packages rcompanion, vcd, psych, desctools, and epitools. Finally, PCA analysis was performed using packages FactoMineR (80) and factoextra (81), taking into consideration 71 taxa at species level.
Data availability. These sequence data have been submitted to the NCBI repository Sequence Read Archive (SRA) databases under accession numbers SAMN16442373 to SAMN16442378; SAMN16442391 to SAMN16442396; SAMN16442397 to SAMN16442402; SAMN16442409 to SAMN16442414; SAMN16442427 to SAMN16442432; SAMN16442444 to SAMN16442449; SAMN16442450 to SAMN16442455, and SAMN16442462 to SAMN16442467, with BioProject number PRJNA669646. Supplemental data, including Exel files of elaborated data obtained from qPCR for target microbial groups and ARGs and NGS data categorized at phyla, family, and genera levels, are available on reasonable request from the corresponding author.