Differences in honey bee bacterial diversity and composition in agricultural and pristine environments – a field study

Agrochemicals and biocides are suspected to cause a dysbiosis of honey bee microbiota, decreasing colonies ability to respond to the environment. As a first step to investigate agriculture and beekeeping impact, hives bacteriomes from an anthropized environment (Agri-env) were compared to that of pristine’s (Prist-env). 16S rRNA sequencing evidenced differences in richness and composition between sample types (Gut (G), Brood (B), Bee-bread (BB)) and environments. Higher opportunist loads and shifts toward taxa capable of metabolizing insecticides were observed in G and B at Agri-env, while beneficial bacteria were enriched in Prist-env. Bacteria in BB did not differ, the acidity of the niche outweighing the influence of external factors. Results showed the environment plays a major role in shaping honey bee microbiota, the agricultural realm inducing a bacterial disruption that would let to colonies vulnerability. In contrast, a less susceptible bee will be promoted in less anthropized locations.


INTRODUCTION
Honey bees (Apis mellifera ) are essential for global agriculture and are key players in many natural ecosystems, being the most important pollinators and honey producers. Unfortunately, sudden colony losses have occurred worldwide, which raised the need to explore the factors behind such dramatic decline (Simon-Delso et al. 2014). There are evidences of several factors being involved, such as the decrease in genetic diversity, habitat loss, poor nutrition, the presence of a wide range of pathogens and Electronic supplementary material The online version of this article (https://doi.org/10.1007/s13592-020-00779-w) contains supplementary material, which is available to authorized users. parasites, the introduction of alien species, climate changes, or intensive agriculture practices (vanEngelsdorp and Meixner 2010;Crotti et al. 2013;Goulson et al. 2015;DeGrandi-Hoffman et al. 2017). The combination of these factors has resulted not only in large losses of colonies but also in an increase of honey bees' susceptibility to parasites (such as to Varroa destructor ). The latter has led to the intensification of treatments against parasites to ensure the survival of the colonies, which, together with the pesticides / insecticides used in agricultural systems to combat pests, are now believed to be the two main causes responsible for high losses of colonies (Kakumanu et al. 2016;Hu et al. 2017).
Multiple pesticide residues have been detected in bees, wax, pollen, and bee-bread. Mullin et al. (2010) found alarming levels of miticides and agricultural pesticides in honey bee colonies from the US and Canada, with around 98 insecticides, fungicides and herbicides detected in pollen. These compounds can alter host reproduction, development, or the immuno-competence of bees, as well as their neurophysiology. Thus, they impact colony health and bees' susceptibility to parasites / pathogens colonization and growth (Thompson 2003;Desneux et al. 2007;Anderson et al. 2011;Staveley et al. 2014). In addition, there is ample evidence of the negative effect of some common pesticides on honey bees gut microbial consortia. For instance, the application of neonicotinoid pesticides and their relationship with pollinators decline have been extensively debated (Suryanarayanan 2015;Mitchell et al. 2017). Neonicotinoids can have synergistic effects with some pathogens, e.g. Nosema (Alaux et al. 2010), increasing its damage on honey bee health. Similarly, Chlorothalonil and Coumaphos have been shown to alter the relative abundance of some beneficial bacterial groups (Lactobacillales and Bifidobacteriales among others) as well as in the variation of the functional profile of some honey bee gut taxa (Kakumanu et al. 2016).
The microbiota associated to honey bees has received great attention in last years, with studies trying to address the symbiotic and pathogenic associations to better understand their intrinsic relationship with bee health (Cox-Foster et al. 2007;Evans and Spivak 2010;Engel et al. 2016). Specifically, some studies have shown that the honey bee gut core is dominated by few phylotypes, 7-12 bacteria that occur regularly, such as strains of Acetobacteriaceae (Alpha 2), Bifidobacterium, Gilliamella (Gamma 1), Frischella (Gamma 2), Snodgrasella (Beta), Enterobacteraceae, Lactobacillus spp., and other Firmicutes Anderson et al. 2013;DeGrandi-Hoffman et al., 2017). Interestingly, honey bee, and social insects in general, have significantly fewer immune genes than expected (Hu et al. 2017) and thus, it has been proposed that their associated microbiota supports important functions. In this regard, several studies have already shown the ability of some honey bee microorganisms to inhibit the progression of diseases caused by bacteria and fungi (Reynaldi et al. 2004;Anderson et al. 2013), to stimulate host innate immune response (Evans and Lopez 2004) or to help in bees nutrition (Crotti et al. 2013), therefore playing an important role in honey bees disease and pathogen / parasite susceptibility (Evans et al. 2006).
Biocides applied in crops and over hives have been shown to alter the diversity and structure of the microbial community associated with the gut and bee-bread of honey bees (Crotti et al. 2013;Simon-delso et al. 2014;Jones et al. 2018), however the majority of the experiments are conducted under a controlled exposition of pesticides or miticides in experimental hives, in contrast to the dosage variability and cumulative exposition occurring in field hives. In addition, biocides transformation within the hive is yet unknown, and their interactions with beneficial microorganisms (Mullin et al. 2010) is still to be resolved, which would be the key to better understand how agriculture and beekeeping practices affect the health of honey bees and their symbionts. In the same way, those chemicals impact upon the microbiota of other in-hive spaces, besides the honey bee gut, are not yet entirely understood. To address this knowledge gap, in the present study, the bacterial community associated with honey bee colonies (gut, bee-bread and brood), their bacteriome, located in two environments of contrasting anthropic pressure have been characterized through 16S rRNA amplicon sequencing.
The bacteriome of a big commercial apiary surrounded by crops and continuously treated against mites has been compared with that of honey bees inhabiting in a pristine, nearly unpopulated island, where colonies survived without any human intervention since 2012 despite the presence of mites. Honey bees in the commercial apiary studied here are affected by more intensive agricultural practices, such as herbicides/pesticides application as part of pest management in crops, monocultures that potentially reduce the quality and quantity of available nutrients for bees, and /or in-hive miticide applications for controlling pathogen colonization as part of beekeeping practices. We hypothesis that such agricultural stressors compromise colonies health as a result of the elimination and/or growth inhibition of particular beneficial bacteria and the disbyosis of their microbial profile. Collectively, this study is a first step toward a better understanding of the interplay between honey bee microbiota and anthropic stressors with implications for future development of strategies to reinforce the health of colonies.

Sampling locations and management of the colonies
Samples were collected from colonies located in two contrasting sites in Croatia: a) An apiary placed in Čeminac, in the Osječko-baranjska County, the northeast of Croatia (45°40'12.88" N, 18°40'40.84" E) ( Figure S1). This apiary is heavily influenced by agriculture and related practices (herbicides/pesticides, monocultures, etc.). The selected apiary is a commercial one running for twenty-five years, but also contains colonies for scientific research. Colonies are treated every year against the varroa mite with the acaricide Checkmite® (Coumaphos). The apiary is located on the edge of an acacia forest (Robinia pseudoacacia ), surrounded by intensive commercial crops, mainly rapeseed, wheat, sunflower, corn and soybeans, as well as fruit trees of apple and plum. This environment, more influenced by agricultural activities, will be referred to as"Agrienv". b) An apiary placed in Unije island (16.92 Km 2 ), in the Adriatic Sea (44°38'59 "N, 14°1 4'44" E) ( Figure S1). The island, with only 90 inhabitants, can be considered a pristine wild place, offering to honey bees an environment completely isolated from practices and products associated with the intensive agriculture or beekeeping management. The minimum distance from the island to the mainland is around 40 Km that ensures honey bees isolation from the contact with intensive crops as they can only travel an average of 4.5 Km during their foraging flights (Seeley 1985). The colonies located in Unije have survived without any beekeeping management since 2012 despite V. destructor and other pathogens being present in these colonies. This environment, hardly influenced by agricultural activities, could be consider a pristine environment and thus hereafter it will be referred to as "Prist-env".

Samples collected
In total, 17 colonies (9 from "Agri-env " and 8 from "Prist-env") were sampled during the spring of 2016 (Table I). Three sample types were collected at each colony: 1) a 50 mL falcon tube filled with worker bees, 2) a 12-15 cm 2 piece of brood in a zip bag, and 3) an 8 cm 2 piece of bee-bread in a zip bag. All the material used during the sampling was sterilized with ethanol 100% and ultraviolet light to ensure sterile conditions. After taking samples from one colony, all the beekeeping tools were sterilized and clean gloves were used for the next sampling. Collected samples were frozen in dry ice in the field and the cold chain (-20°C) was kept until their arrival to the Genetics laboratory of the University of the Basque Country, Spain, where they were stored at -80°C until DNA extraction.

DNA extraction
For gut samples (G), 10 bee guts were extracted by dissection. Guts were placed in a 1.5 mL tube and 600 μl of PBS 1x were added to strongly homogenize it using a Precellys24 tissue homogenizer (Bertin, Montigny-le-Bretonneux, France). Then, it was centrifuged for 1 min at 13,000 rpm and the supernatant ("clarifiat") was placed in other 1.5 mL tube. This process was repeated one more time by adding 400 μl PBS 1x in order to recover the maximum sample possible without taking solid remains. The resultant supernatant was added to the previously recovered and the total was stored at -20°C. 200 μl of "clarifiat" was used as input for DNA extraction with the QIAamp® DNA Mini Kit (Qiagen) following the indications specified by the manufacturer.
Regarding brood samples (B), three black eyes pupae were pooled per colony, removing the head (due to the known PCR inhibitors (Boncristiani et al. 2011)). In addition, the brood cells of the used pupae were cleaned with PBS 1x and this liquid was added to the brood "clarifiat". From this "clarifiat" 200 μl were taken to start the DNA extraction following the QIAamp® DNA Mini Kit protocol.
Regarding bee-bread samples (BB), three bee-bread cells were randomly chosen and their content (bee-bread balls) was placed in a 1.5 mL tube. 200 μl PBS 1x, 400 μl of AL buffer, and 40 μl of Proteinase K was added. Cell lysis was performed by incubating the mixture at 56°C during 1 hour with 600 rpm, shaking it in a vortex each 15 min. After the incubation, 600 μl of Phenol-Chloroform-Isoamil alcohol (25:24:1) was added to the liquid, it was shaken and centrifuged for 15 min at 14,000 rpm. The supernatant was recovered and placed in a new 1.5 mL tube and 600 μl of Chloroform was added, then shaken and centrifuged during 5 min at 14,000 rpm. The supernatant was recovered and placed in a new 1.5 mL tube, 400 μl of ethanol 100% were added into it, and the tube was shaken. All the liquid was then placed in the columns provided in QIAamp® DNA Mini Kit. The extraction was carried on from this step following the manufacturer's instructions.

PCR amplification and bacteria sequencing
To determine the bacterial community present in the samples, the amplification of the V4 region of the 16S rRNA gene, located in the SSU, was performed using primers 515F-806R ("Earth Microbiota Project"; http://www.earthmicrobiome. org/) that contained a barcode sequence (12 bp) in the forward primers. The PCR conditions for the V4 region amplification were the following: 0.3 μL Taq polymerase (Go Flexi Promega Taq), 5 μL of Taq Buffer (5x), 2 μL MgCl 2+ (50mM), 2 μL of 2.5 mM dNTPs, 0.5 μL each primer (10 μM), 1 μL extracted DNA (15-60 ng/μL), and water up to 25 μL of the total volume. The PCR cycles used were: an initial denaturing step at 95°C for 4 min; then 35 cycles of 15 sec at 95°C, 30 sec at 50°C, 30 sec at 72°C; and a final elongation during 2 min at 72°C. The above explained conditions were used for the amplification of the three types of samples (G, B, and BB). For BB this later protocol will be referred to as "without blocking primers" test.
Additionally, in BB samples, a second amplification protocol was tested. As the primers used in the present study (515f-806r) also amplify chloroplast and mitochondria, a PCR including PNA clamps (Lundberg et al. 2013), called "blocking primers" in this article, was tested to block the chloroplast and mitochondrial DNA amplification potentially coming from flower's pollen source. PCR reagents were added in the same concentrations as previously described, but in this test 0.35 μL of the PNA clamps (mPNA and pPNA in 1:1 Table I. Colonies sampled during spring 2016 in the two locations studied: Unije Island (Pristine environment) and Cěminac apiary (Agriculture environment) proportion (50 μM each one)) was added to the PCR mix. In this case the PCR cycles were an initial denaturing step at 95°C for 45 sec; and 34 cycles of 15 sec at 95°C, 10 sec at 78°C, 30 sec at 50°C, 30 sec at 72°C. This protocol will be referred as "with blocking primers" test in the rest of the article. All PCR products were checked on a 1.5% agarose gel stained with ethidium bromide.
DNA purification of the three types of samples was performed with the CleanPCR kit (Cleanna), using magnetic beads. Samples were quantified by Qubit TM v2.0 (ThermoFisher Scientific) and normalized in 8 pM pool. Paired-end sequencing of the pools was carried out on an Illumina MiSeq sequencer at the Sequencing and Genotyping Unit of the University of the Basque country (SGIKER), with the kit v2 PE 2 x 150 bp (300 cycles), adding 10% of PhiX.

Quality control, processing and taxonomic assignment of the 16S rRNA gene
The qualities of the forward and reverse raw sequences were checked with Sickle v1.33 (Joshi and Fass 2011) using a threshold of Q20, and the assembly of forward and reverse sequences was performed with Pear v0.9.10 program (Zhang et al. 2013), using an overlap of 15 bp. The fastq-barcode.pl script (Smith 2012) was used to remove non-existent barcodes from the fastq achieved in the previous step by Pear. Then, sequences taxonomic assignment was carried out using QIIME v1.9 (Caporaso et al. 2010a), through the script pick_open_reference_otus.py against Silva 128 database version (https://www. arb-silva.de) at 97% similarity level. Briefly, this script wraps UCLUST (Edgar 2010) for clustering the sequences in Operational Taxonomic Units (OTUs), with 97% of similarity among the sequences (in the present study), PYNAST (Caporaso et al. 2010b) to align OTU sequences, and FastTree (Price et al. 2010) to construct a tree with the sequences. Mitochondrial and chloroplast sequences were removed using the script filter_ taxa_from_otu_table.py of QIIME v1.9. All OTUs with less than 10 sequences were also removed, using the script filter_otus_from_otu_ table.py . Finally, the dataset was split by sample type.
After all the filtering's, the dataset for gut samples contained a total of 2,598,419 reads (mean/colony = 99,939), 305,310 reads were found in brood samples (mean/colony =33,923), while the dataset for bee-bread contained 393,707 reads (mean/colony = 19,685), and increase to 630,848 reads (mean/colony = 31,542) when blocking primers were used for library construction.
2.6. Characterization of the microbial diversity and composition per sample type: gut, brood and bee-bread The diversity present within the three sample types was calculated to obtain alpha and beta diversity estimates using R (phyloseq package). Previously the dataset was normalized using normalize_table.py script and DESeq2 method in QIIME1.9, that is the recommended analysis when fewer than 50 samples are to be studied. The taxa present in each sample type and their abundance were collapsed by their mean, through QIIME script collapse_samples.py and the results were represented in the bar plot by class taxonomy level in order to summarize the data. An ANOSIM test was calculated comparing the beta diversity among the different sample types using the QIIME script compare_categories.py , to test whether there are significant differences among the bacterial communities associated with gut, bee-bread and brood sample types.

Assessing microbial differences between environments: Agri-env vs Prist-env
The analysis was performed separately for each sample type. The gut dataset included 17 samples, the brood dataset was composed of 9 samples and 11 colonies were used in the BB analysis (Table I). Alpha rarefaction and alpha diversity indexes were calculated with QIIME v1.9 and R v3.3.3 phyloseq package respectively. In order to test whether the microbial diversity differed between anthropic and pristine environments, alpha diversity indexes were estimated grouping samples by environment and whenever it was possible the Mann-Whitney-Wilcoxon test was performed to get the significance of the comparisons. Beta diversity was calculated using the B r a y -C u r t i s d i s t a n c e s m e t r i c s w i t h beta_diversity.py script (QIIME) and this matrix was used to perform an UPGMA tree (performed with QIIME), and a Principal Coordinate analysis (PCoA) plot (performed with R phyloseq package) for visualizing the similarities/dissimilarities of the microbiota structure in the samples between the two environments. The statistical significance of grouping samples by environment was checked by ANOSIM test.
Kruskal-Wallis test was carried out with the script group_significance.py (QIIME) to identify bacterial groups that were more relevant in the microbiota structure differences between the two environments. The abundances of the relevant bacterial groups, at family and genus level, were studied and represented in bars plots through the QIIME script summarize_taxa_through_plots.py . A heatmap using pheatmap package of R was constructed with the OTUs (with higher abundances than 0.1%) that where significant in the Kruskal-Wallis for "environment" category.

Assessing blocking primers efficiency in BB samples amplification
Since high chloroplast and mitochondria presence was expected in bee-bread, the amplification and sequencing of those samples were performed using two protocols: "without blocking primers" and "with blocking primers" (including PNA Clamps in the library construction step). In order to obtain more robust results when comparing both tests, 9 extra samples of BB were collected in the same apiaries and colonies during the next fall (2016) and added to the analysis, making a total of N=20 samples.
The abundance of chloroplast sequences was estimated for each sample and protocol. The community structure and composition obtained with the two protocols was also compared. Briefly, those analyses included alpha rarefaction and alpha diversity estimates, normalization, and ordination and clustering representations, together with significance tests using QIIME scripts and R packages (as defined in the previous section). The alpha diversity indexes were estimated by grouping samples by protocol, and a Mann-Whitney-Wilcoxon test was performed for evaluating if the alpha diversity differences between the protocols used were significant. Families abundances per protocol was visualized in bar plots, and beta diversity matrix was calculated and used to perform a PCoA plot to check whether samples clustered according to the amplification protocol used. ANOSIM test was carried out to address if the difference between the bacterial communities recovered by each method was statistically significant. Finally, Kruskal-Wallis test was performed to identify if any genus detected with the protocol "without blocking primers" was missing from the ones recovered in the protocol "with blocking primers", and thus, failed to be amplified when adding the blocking primers in the PCR.

Characterization of bee gut, brood and bee-bread microbiota
Bacterial richness as estimated by Chao-1 index differed among the three sample types ( Figure 1a). Overall, gut samples (G) had the lowest but most homogenous alpha diversity values. On the contrary, Bee-Bread (BB) showed the highest alpha diversity values with the highest heterogeneity in richness values. Beta diversity plot ordination showed clear community composition differences among the three sample types (Figure 1c) and ANOSIM test for "sample type" variable was significant (P = 0.001; R = 0.992).
The majority of the taxa detected in G belonged to Firmicutes phylum and four bacterial classes, Gammaproteobacteria, Betaproteobacteria, Alphaproteobacteria, and Actinobacteria (Figure 1b), showing a very stable bacterial community. Brood samples (B), however, showed several taxa in very low proportion and numerous taxa appearing in only one or two samples (see Figure S2), resulting in a high heterogeneity in both, richness and composition (Figure 1a and c). The dominant classes in B were Alpha-, Gamma-, and Delta-proteobacteria (at 14.5%, 13%, and 10.6% abundance, respectively). Flavobacteria or Betaproteobacteria reached abundances of 5.3% and 4.6% (Figure 1b). The rest of bacterial groups showed minor abundances (< 3% for all sample types). Similar to brood samples, BB were also composed of high number of low abundant taxa (almost all with relative abundances lower than 3%) and few moderately abundant organisms, being Gamma-proteobacteria the most abundant class (14.9%) followed by Alpha-, and Deltaproteobacteria (12.9% and 8.7%, respectively). The majority of families were present in very low proportions (≤ 0.01%) (see Figure S3), although families such as Flavobacteriaceae, Sphingomonadaceae, A n a e r o l i n e a c e a e , R h o d o b a c t e r a c e a e , Planctomycetaceae, Desulfobacteraceae, Alteromonadaceae, Saprospiraceae, and Desulfobubaceae were present in ranges from 5.8% to 1.5%.
The primers used for library construction (515f-806R) are known to amplify chloroplast sequences, and their abundances were as high as 88.5% in some BB samples. In an attempt to avoid chloroplast sequences overwhelming the sequencing run, a modified protocol was tested, that included PNA Clamps in the library construction step ("with blocking primers" protocol). The suitability of the modified protocol was evaluated. Chloroplast detection was considerably reduced to ranges 3.5%-26.8% (more details in Table S1). Alpha diversity estimates were significantly higher ( Fig. 2a; Mann-Whitney-Wilcoxon test (W = 400, P = 6.786e-08 for observed-OTUs values; W = 400, P = 1.451e-11 for Chao-1 index)) detecting higher numbers of different families when PNA Clamps were added (Figure 2b). In fact, the bacterial community was significantly different between protocols according to ANOSIM test (P = 0.001; R = 0.96) and the clustering in the PCoA plot (Figure 2c). The abundance of severa l ge nera (su ch as Methylophilus , Nitrosomonas , Nitrosococcus , Sulfitobacter , Desulfobacterium , Desulfococcus , etc) differed between the two protocols (Kruskal-Wallis test, Bonferroni P <0.05). All of them were detected in lower abundances or not detected in the protocol that did not include PNA Clamps. Thus, for the microbial diversity and composition analysis in BB, the dataset obtained in the modified protocol, referred to as "with blocking primers" protocol, was used.

Gut samples
Honey bee gut samples from the pristine location (Prist-env) had significantly lower Chao-1 values than those in Agri-env samples (Figure 3a, Mann-Whitney-Wilcoxon test W = 9, P = 0.008). However, both environments displayed a similar variance.
Samples were grouped by location in the UPGMA clustering and in the PCoA ordination plot, where two clear groups were observed (Figure 4a). The "environment" variable explained 20.4% of the differences observed in the bacterial community composition (Figure 4b), and those differences were Kruskal-Wallis test identified several taxa belonging to 14 bacterial families with a significantly abundance difference (Figure 5a; details per sample in Figure S4). Lactobacillaceae, Acetobacteriaceae, Neisseriaceae, Sphingomonadaceae, and Orbaceae had high discriminative power to distinguish pristine honey bee gut samples, while higher abundances of Enterobacteriaceae (particularly Pantoea and Enterobacter genus representing the biggest differences) and Bifidobacteriaceae were characteristic in the agricultural location (Agri-env).
Within those 14 families (Figure 5b; details per s a m p l e i n F i g u r e S5 ) , L a c t o b a c i l l u s , Commensalibacter , and Snodgrassella genus presented high differences between the two environments, being particularly abundant in the Pristenv gut samples. Others, such as Sphingomonas , Enterococcus , Gilliamella , Methylobacterium, Frischella , Saccharibacter , Acinetobacter , and Pseudomonas , showed slighter abundance differences but were also significantly enriched in colonies from the non-anthropized condition.
A clear division among pristine and agricultural samples was evidenced in a heatmap constructed with the 118 OTUs included in the differentially abundant genera detected through Kruskal-Wallis ( Figure 6, OTU clusters details in Supplementary Table S2). Clusters formed by one (25 and 26) or more than one (2, 9, 13, 16, and 22) specific strains of Enterobacteriaceae genera were highly abundant in Agri-env, being their abundances very low in Prist-env. Bifidobacterium genus (8 and 10) cluster followed the same trend. On the contrary, strains of Lactobacillus (cluster 1, 5, 12, 24) or Commensalibacter (cluster 4 and 30) were mostly present in pristine samples. Particular strains belonging to Acinetobacter , Snodgrasella , Gilliamella , Frischella , or Sphingomonas , also considered beneficial in the honey bee gut microbiome, appeared in both environments. However, specific OTUs/strains within them appeared to be more associated with one or the other location (e.g.: clusters 3, 11, 15, 21, and 27 in Agri-env samples; clusters 7, 20, and 28 in Prist-env).

Brood samples
Prist-env colonies showed higher alpha diversity values and higher variability than agri-env colonies (Figure 3b). When accounting for OTUs present in at least 50% of the samples, all the analysis performed showed a clear differentiation between locations: 1) The UPGMA tree grouped both types of samples at different branches (Figure 7a). 2) The PCoA plot showed a clear grouping for the Agri-env samples and Axis 1 explained 48.6% of the variability (Figure 7b), and 3) ANOSIM test for "environment" resulted in a significant (P = 0.012) and high R value (R = 1) after M45 colony was removed (a colony shown to have a very dissimilar microbial composition).
16 bacterial groups and 4 archaea were found to have significantly different abundances at a family taxonomic rank (Figure 8, details for each sample in Figure S6). The biggest differences in abundance were found for members of Sphingomonadaceae and Methylobacteriaceae families (Sphingomonas and Methylobacterium genera accounting for highest differences (Figure 9)), all more abundant in Agri-env. Various members within Flavobacteriaceae family, specifically, Flavobacterium and Chryseobacterium were also overrepresented. Prist-env showed higher abundances of Anaerolineaceae, Alteromonadaceae and Archaea families (Figure 8), as well as an enrichment of Ignavibacterium , Haliea , Halioglobus , and Pseudohaliea genera ( Figure 9). Actibacter , some unassigned Flavobacteria (Figure 9), and two Robginitalea genus OTUs were also more abundant in the pristine non-human manipulated colonies.

Bee-bread samples
Alpha diversity values were very variable within the colonies from both environments, but particularly in Prist-env, similar to the results found in brood samples (Figure 3c). The bacterial community did not demonstrate a grouping related to the environment, appearing mixed in the UPGMA and PCoA plot (Figure 10a and b; non-significant ANOSIM value).

Characterization of bee gut, brood and bee-bread microbiota
It is already known that the diversity of the honey bee gut is modest and very stable, with a core microbiota shared across individuals and studies . The results obtained in the present study were in line with this evidence: gut samples had the most homogeneous alpha diversity values across all sample types studied (Figure 1a), and the bacterial composition found was in concordance with those identified in previous studies. Gammaproteobacterias and Firmicutes were the most abundant (Figure 1b), with Gilliamella (6.1%), Frischella (4.6%) and Lactobacillus (28.5%) being the main genera within those groups. Betaproteobacteria were also abundant, being Snodgrasella genus (12.3%) the  one with the highest representation. The repeated presence of these bacterial groups across several studies that account for different honey bees populations of diverse genomic background, collected at different seasons and regions, lead hypothesize their implication in central functions, being key for honey bee and colony health, as it has been suggested for other animals or humans (Kau et al. 2011;Martinson et al. 2012;Moran et al. 2012).
Brood and Bee-bread microbiota have been yet scarcely investigated. Both samples are intrinsically related to each other, and partially related to the bee gut. Larvae are fed by nurses with glandular secretion, pollen bread and honey (Winston 1985), being this their first contact with the microbiota present in adult worker guts and beebread. After capping, there is no more contact from the exterior. Therefore, the pupae's microbiota as sampled in this study, should in principle resemble the larvae microbiota, that would be a mixture of bee-bread, honey and worker gut microbiota, as well as the cell environment. On the other hand, Bee-Bread is a mixture of flower pollen, nectar, salivary of the honey bee hypopharyngeal glands, and the comb/cell environment (Anderson et al. 2014). The bee-bread is, thus, expected to reflect the microbiota coming from the surrounding environment (e.g. crops) and from bee saliva and gut. The results from the present study showed that BB and B had most of their bacterial groups in common (e.g. Gamma-, Delta-, or Alphaproteobacteria, Flavobacteria or Anaerolineae) and that this niches are composed of high abundances of groups thought to be noncore or sporadic bacteria in gut, such as Desulfobacterales (Deltaproteobacteria) and Flavobacteriaceae (as previously identified by Newton andRoeselers (2012), or Olivieri et al. (2012)), but their functions in the hivemicroenvironment is unknown yet.
Despite the similarities between B and BB, some groups showed higher abundances in brood samples, such as Comamonadaceae (3.4% vs. 1.2%), Desulfobacteraceae (3.9% vs. 2.9%), Methylobacteriaceae (1.7% vs. 0.6%), Desulfobacterales (5.8% vs. 4.5%), or the Actinobacteria Corynebacterium (0.8% vs. 0.3%). The latter is known to inhibit fungal growth (Anderson et al. 2013), and thus, it is plausible that it may help to protect bee brood against fungal pathogens. This would support the hypothesis that the scarcity of bacteria in larvae (and early developmental stages within brood) could be due to the brood cell containing growthsuppressing antimicrobial agents (Martinson et al. 2012).
The microbial profile found for BB showed a specialization towards acidic environments with the enrichment of acid-tolerant microorganisms (e.g. Acidobacteria, Actinobacteria, some Flavobacteria, etc). This is in accordance to Anderson et al. (2014), that suggested that there is an initial inoculation of microbiota from bees during collection and consolidation of corbicular pellets, but that afterwards, due to the acidity of the environment and high content of sugars, the hive stored pollen bacteria are selected towards microorganisms adapted to such condition.
The abundances of Acetobacteraceae (0.5%), Actinobacteria (4.2%), Enterobacteriaceae (1.2%), Actinomycetales, Fructobacillus or Staphylococcus (< 1%) in BB samples from the present study were lower than expected in comparison to what was found by other authors (Anderson et al. 2013(Anderson et al. , 2014. These differences could be attributed to various factors associated with the season and/or region investigated, that would determine the quality and type of the pollen available for bees. Additionally, methodological differences across studies could also be behind some of the differences. In fact, in the present study we found differences in the richness and Figure 6. Heatmap clustering of the bee gut OTUs with significant abundance differences between the two environments. OTUs with relative abundances higher than 0.1% and belonging to the 20 genera that showed a significant result in the Kruskal-Wallis test are represented. The heatmap was constructed using a Manhattan matrix. Further details of the members comprising each cluster are available in Table S2. taxonomic profiles of the same sample sequenced by using two different library construction protocols; with and without including PNA clamps (Lundberg et al. 2013) to block the amplification of chloroplast and mitochondrial DNA potentially coming from pollen. High reduction of the unspecific non-target chloroplast DNA amplification was observed with the addition of PNA Clamps. Importantly, this modified protocol increased the detection power of the families present in BB samples. Moreover, higher alpha diversity values were recovered, and any taxonomic group amplified in the protocol "without blocking primer" was also detected when blocking primers were used. All the above results demonstrated that the addition of PNA Clamps constitutes a methodological improvement to investigate plant source samples associated with bees, in particular when  Relative abundances (mean percentage) and standard deviation (error bars) of A) bacterial families (kruskal wallis p-value < 0.05), B) genera showing a differential distribution within the significant families. Further details for each sample are available in Figure S6. using primers that do not discriminate between plastids and bacterial organisms.

Gut
It is known that the environment, the diet, age, etc., modifies the gut microbiota (Martinson et al. 2012), as well as several mechanisms of bee behavior, such as trophallaxis, oral-fecal-route, or grooming, that help maintaining it stable across the individuals in the hive (Alberoni et al. 2016). The exposure of bees to the environment around the hive can bring new bacteria, as well as agrochemical substances into the hive. Such compounds can get accumulated and alter bee associated microbiota (Henry et al. 2012;Crotti et al. 2013). Likewise, the contamination of wax with drugs applied in the hives (against varroa mites and other pathogens) have been shown to have a negative effect on microbiota (Tihelka 2018). The derived diversity changes and compositional modifications could leave bees in a vulnerable state with respect to pathogens, since the dysbiosis in the gut could alter the expression of developmental genes and immune system functions, which will decrease the ability of bees to respond to environmental or pathogenic stressors (Kakumanu et al. 2016;Raymann and Moran 2018).
Under this scenario, opportunistic bacteria could take advantage of the microbial shifts, colonizing and proliferating in that niche leading to an increase in richness. This explanation could fit with the greater richness estimates found in the honey bees gut raised in the Agri-env in the present study (Fig. 3a). In this line, an increase of the opportunistic pathogen Serratia marescens was reported by Motta et al. (2018) after the exposure of colonies to glyphosate, which reduced the abundance of symbionts in the community. Moreover, a balanced microbiota is an important component of colony health (Anderson et al. 2011). Therefore, we hypothesized that the differences in microbial abundances found in the present study could lead to differential pathogen susceptibility between the honey bees raised in the two environments investigated. For instance, Pettis et al. (2012) and Wu et al. (2012) observed that the honey bees fed with pollen containing Chlorothalonil were three times more susceptible to Nosema infection than controls. Coumaphos and Chlorothalonil exposure were shown to increase Bifidobacterales and Enterobacteriaceae (Kakumanu et al. 2016). In addition, a decrease in Lactobacillales, considered to inhibit the growth of major honey bee pathogens, was also noticed after the application of Chlorothalonil (Promnuan et al. 2009;Forsgren et al. 2010;Audisio et al. 2011). Likewise, Raymann et al. (2017), showed a reduction in Lactobacillus and Snodgrasella bacteria, among others, as well as an increase of infection by ubiquitous opportunistic pathogens, when antibiotics were applied to honey bees. In concordance, in this study, Lactobacillus and Snodgrasella abundances were depleted and Bifidobacterales and Enterobacteriaceae were enriched in the guts subjected to higher agricultural and beekeeping pressure. Considering that Bifidobacterium is known as a beneficial genus within the gut (Alberoni et al. 2016), its higher abundance under agricultural stressors might have been associated with the decline of other beneficial species that under other circumstances compete in the gut ecological niche with Bifidobacterium. This disruption could alter honey bee metabolisim and/or let the disproportionate growing of opportunistic pathogens like some Enterobacteriaceae . On the contrary, honey bee guts from Unije island had higher proportions of beneficial bacteria, such as Snodgrasella , Gilliamella , Lactobacillus , or Commensalibacter genera. Those bacterial groups were previously proposed to be involved in the protection against parasitic protozoans in bumble bees (Koch and Schmid-Hempel 2011), and thus, they may perform a similar function in honey bees. They also participate in biofilm formation, e.g. Snodgrasella (Betaproteobacteria, Neisseriaceae), promoting fastening and colonization of the gut by other phylotypes such Gilliamella (Gamma-1, Orbaceae) (Martinson et al. 2012). Similarly, Kwong et al. (2017) determined that Snodgrasella alvi and Gilliamella apicola participate in immune gene expression and survival rate after E. coli infection. Frischella , involved in the immune gene expression and melanization response (Emery et al. 2017), showed a weak difference between environments in the present study, but still, were slightly more abundant in Prist-env gut samples. The higher abundances of such beneficial groups Figure 10. Bee-bread samples community dissimilarity based on the Bray Curtis beta diversity matrix. A) UPGMA tree; the colonies from agricultural environment are marked with a black point, B) PCoA with the sample points colored by "environment" and confidence ellipses at 95%. may contribute to the impressively high survival rates of the colonies investigated in the Prist-env, who lasted over 8 years on the island without beekeeping intervention but in the presence of Varroa destructor .
Interestingly, Jones et al. (2018) identified that some Acetobacteriaceae and Lactobacillaceae had higher abundances in honey bee guts raised distant to crops (≥1.25 km away from hives). Similarly, higher abundances of those taxa were found at prist-env samples in the present study. However, the here reported abundances are higher than those found by Jones's and collaborators (2018), probably due to honey bees raised in Unije being located further away from any agricultural system precluding them to be in contact with such stressors. Although we cannot discard that part of the differences observed in the microbiota are due to the island effect and the vegetation present in this place, our results, in combination with Jones's and collaborators (2018), point to the hypothesis that the microbiota associated with the gut of bees is increasingly modified with the intensity of, and distance to, the agricultural environment.

Brood
Brood samples from colonies located in the pristine environment showed higher alpha diversity estimates and higher variability in their bacterial composition, on the contrary to the results found for honey bee gut. Interestingly, microbiota composition differences were found according to the environment, that would be in line with the assumption that brood is mainly composed by environmental phylotypes, and therefore, dependent on the environment within and outside the hive.
Prior to the pupae stage and cell capping, larvae are fed by nurse workers with glandular secretion, pollen and honey (Winston 1985). Therefore microorganisms that are transferred to the brood during this process would remain within the cell. In such scenario, potential modifications on the herbage, the bacterial composition in pollen surface, or bee gut caused by within-hive miticides and/or pesticides are expected to be reflected in the bacterial community present in brood cells.
Accordingly, members of the Sphingomonas family, which have enzymes able to degrade some insecticides components, such as organochlorines, or synthetic pyrethroids (Russell et al. 2011), were higher in samples from the commercial apiary. Conversely, the environmental bacteria A n a e r o l i n e a c e a e , A l t e r o m o n a d a c e a e , S a n d a r i n a c e a e , I g n a v i b a c t e r i a c e a e , Acidimicrobiaceae, or Archaea were clearly more abundant in the pristine environment. To what extent this microbial differences are related to honey bee colony health is still to be resolved.
A differential abundance of the varroa mite could have also in part enlarged the differences among the colonies and between the environments. Honey bees in Unije did not receive any acaricide treatment but were infected with varroa mites. In this regard, Flavobacteria showed abundance differences between the environments. This bacteria has been previously identified as part of ticks and varroa mite microbiota by Sandionigi et al. (2015), thus, Flavobacteria might play a role in the final pathogenic activity of the parasitization and therefore might be involved in the differential success of varroa parasites in those colonies. Further studies focusing on the microbial changes along different stages of larvae and pupae, together with the study of the brood cell environment, are necessary to shed light into the function of these and other bacteria during disease infection and parasitization.

Bee-bread
Pollen pellets collected by bees are packed into wax storage cells within the hive, and become bee-bread through fermentation. BB is a product rich in microbes, vitamins, lipids and proteins (Anderson et al. 2013). As in the case of brood, BB samples from the pristine environment showed greater richness and variability. The reduced diversity in the commercial apiary might be explained by the agricultural environment having a more homogeneous vegetation, usually limited to a few crop species and limited availability of wild plants. In addition, as herbicides and pesticides are applied routinely into these systems, their application could eliminate and/or constrain the growth of certain microorganisms (Aktar et al. 2009). Furthermore, while agrochemicals keep in control insects and parasites, they also alter plants epiphytic community (Mullin et al. 2010;Klein et al. 2017). In this sense, Anderson et al. (2013) suggested that changes in the pollination environment, natural or anthropic, could alter the evolution, abundance, transmission rate and/or survival of microbes present in the floral nectar, phyllosphere and local water sources. Those microbial shifts might ultimately impact the hive microbiota, as bees will act as vectors of plant microbiota into the hive environment during pollen and nectar foraging (Graystock et al. 2015).
However, despite the higher richness showed by the pristine Bee-Bread samples, and unlike to the findings in Gut and Brood samples, the bacterial community composition in BB did not differed between the environments investigated. Several reasons could be behind the absence of such differenciation. On the one hand, the pressure exerted by the acidity of the bee-bread media might have outweighted the potential influence of external factors, such as pesticides accumulation, the landscape, the season, etc. On the other hand, it could have a technical explanation. In the present study each BB sample was composed by microbial DNA coming from just three bee-bread cells, that were randomly chosen per colony. Considering the high within-sample variability found for BB sample type, it is plausible that the small number of comb cells studied per colony might have limited our ability to obtain a representative diversity of the whole hive. Therefore, it cannot be discarded that the possible microbial differences between the environments might have gone hidden in bee-bread samples, and thus, further studies targeting a higher number of cells within each colony, in combination with a palynological analysis to clarify the effect of plant composition, is recommended.

CONCLUSIONS
The honey bee associated microbiota is believed to play a key role in bee health regulation (Jones et al. 2018). The application of agrochemical substances or in-hive miticides associated to agricultural environments (intensive crop farming and beekeeping management activities) can result in microbial dysbiosis with diversity and composition modifications that may increase bees' vulnerability to other abiotic and biotic stressors. In the present study, the microbiota of honey bee gut, brood and bee-bread (food storage) has been investigated comparing the microbial profiles of colonies raised under constant agricultural stressors (Agri-env) with that of a pristine insular environment where colonies survived without any human intervention despite the presence of mites (Pristenv). We cannot establish causality between differences in the microbiota composition with pesticide/miticides exposure and colony susceptibility to pathogen infection, since the landscape or the island effect also contribute to the bacteriome differences found. However some of the results obtained head in such a direction and could be partially attributed to the agricultural pressure. For instance, a shift in bacteria capable of metabolizing insecticide compounds and in taxonomic groups previously found to respond to commonly used pesticides was observed in honey bee gut and brood samples located in the environment with higher agricultural and beekeeping pressure (Agri-env). In addition, an enrichment of microorganisms shown to be beneficial, due to their role in immune gene expression and/or involved in the protection against parasites, were found in Prist-env gut and brood samples. Interestingly, those colonies showed high survival rates even without beekeeping management pointing to the pristine environment leading to a healthierstronger bee. This study is a first promising step to resolve the role of the microbiota within the hive and the impact of abiotic factors, which deserves to be further investigated. More studies including long-term experiments and a functional approach, rather than just based on taxonomic characterization, are still needed for a deeper understanding of the bacterial role in the colony health.