Microbial Community Structure among Honey Samples of Different Pollen Origin

Honey’s antibacterial activity has been recently linked to the inhibitory effects of honey microbiota against a range of foodborne and human pathogens. In the current study, the microbial community structure of honey samples exerting pronounced antimicrobial activity was examined. The honey samples were obtained from different geographical locations in Greece and had diverse pollen origin (fir, cotton, fir–oak, and Arbutus unedo honeys). Identification of honey microbiota was performed by high-throughput amplicon sequencing analysis, detecting 335 distinct taxa in the analyzed samples. Regarding ecological indices, the fir and cotton honeys possessed greater diversity than the fir–oak and Arbutus unedo ones. Lactobacillus kunkeei (basionym of Apilactobacillus kun-keei) was the predominant taxon in the fir honey examined. Lactobacillus spp. appeared to be favored in honey from fir-originated pollen and nectar since lactobacilli were more pronounced in fir compared to fir–oak honey. Pseudomonas, Streptococcus, Lysobacter and Meiothermus were the predominant taxa in cotton honey, whereas Lonsdalea, the causing agent of acute oak decline, and Zymobacter, an osmotolerant facultative anaerobic fermenter, were the dominant taxa in fir–oak honey. Moreover, methylotrophic bacteria represented 1.3–3% of the total relative abundance, independently of the geographical and pollen origin, indicating that methylotrophy plays an important role in honeybee ecology and functionality. A total of 14 taxa were identified in all examined honey samples, including bacilli/anoxybacilli, paracocci, lysobacters, pseudomonads, and sphingomonads. It is concluded that microbial constituents of the honey samples examined were native gut microbiota of melliferous bees and microbiota of their flowering plants, including both beneficial bacteria, such as potential probiotic strains, and animal and plant pathogens, e.g., Staphylococcus spp. and Lonsdalea spp. Further experimentation will elucidate aspects of potential application of microbial bioindicators in identifying the authenticity of honey and honeybee-derived products.


Introduction
The honey market constitutes an important agricultural sector, with the global honey production being estimated at approximately 6.6 billion USD in 2015 [1]. Moreover, honeybees are essential for agriculture, ecology, and the maintenance of life, as the pollination assures the reproduction of plants; thus, beekeeping is actually an issue of concern all around the world [2]. In the European Union (EU), the second largest producer of honey Thus, the aim of the present study was to identify, for the first time, the predominant microbial communities in Apis mellifera honey samples of various pollen origin (fir, cotton, fir-oak, and arbutus) exhibiting high antimicrobial activity against common foodborne and human pathogens, and to comparatively evaluate their microbial community structure and ecological indices, through the application of high-throughput sequencing techniques.

Results and Discussion
Four honey samples from bees fed with pollen and nectar of various melliferous plant species were examined in terms of their microbial community structure. Regarding diversity indices, the fir and cotton honeys exhibited significantly greater diversity than the firoak and Arbutus unedo honeys (a < 0.05, in Duncan's multiple tests for Chao1, Fisher, Shannon, and Simpson diversity indices) (Figure 1). No statistically significant differences were identified between fir and cotton honey, or between fir-oak and Arbutus unedo honey regarding all ecological indices estimated, except for the Shannon index, where Arbutus unedo honey showed the least score ( Figure 1).

Figure 1.
Diversity indices of the examined honey samples. Bars represent ± standard errors of means. Analysis of variance (ANOVA) using Duncan's multiple post hoc tests at a significance level of 5% (a < 0.05) were carried out to identify statistically significant differences. No letter in common at the top of the subgraphs is indicative of statistically significant differences.

Figure 1.
Diversity indices of the examined honey samples. Bars represent ± standard errors of means. Analysis of variance (ANOVA) using Duncan's multiple post hoc tests at a significance level of 5% (a < 0.05) were carried out to identify statistically significant differences. No letter in common at the top of the subgraphs is indicative of statistically significant differences.
The examined fir honey was dominated by lactic acid bacteria of the genus Lactobacillus (19.80 ± 0.59% of the total relative abundance), followed by members of the genera Bradyrhizobium (11.93 ± 0.20%) and Pseudomonas (11.42 ± 1.68%) ( Figure 2). Specifically, L. kunkeei (current taxonomic name Apilactobacillus kun-keei) represented 19.34 ± 0.54% of the total relative abundance at species level in fir honey (Table 1).   Lactobacilli appeared to be favored in honey produced from fir-originated pollen and nectar, with Lactobacillus population being more pronounced in fir honey compared to fir-oak honey, indicating a proliferation of this taxon under feeding of honeybees with fir ( Figure 3). Interestingly, Bradyrhizobium has been recently reported to exert antagonistic activity against various pathogens [23]. Recently, co-existence of lactic acid bacteria and Bradyrhizobium resulted in enhanced nodulation with positive impact on plant growth, a fact that may be indicative of co-evolution of these microbiota in certain honeybee hosts [24].  Lactobacilli appeared to be favored in honey produced from fir-originated pollen and nectar, with Lactobacillus population being more pronounced in fir honey compared to firoak honey, indicating a proliferation of this taxon under feeding of honeybees with fir ( Figure 3). Interestingly, Bradyrhizobium has been recently reported to exert antagonistic activity against various pathogens [23]. Recently, co-existence of lactic acid bacteria and Bradyrhizobium resulted in enhanced nodulation with positive impact on plant growth, a fact that may be indicative of co-evolution of these microbiota in certain honeybee hosts [24]. The predominant taxa in cotton honey were Pseudomonas, Streptococcus, Lysobacter, and Meiothermus, representing 14.35% ± 2.10%, 7.66% ± 1.69%, 7.07% ± 1.75%, and 6.53% ± 3.93% of the total relative abundance, respectively ( Figure 2). Pseudomonas is an inhabitant of bee pollen, contributing to the decomposition of pollen walls [25]. As a result, this taxon is part of honeybee microbiome, which is commonly detected in honey [25][26][27]. Notably, streptococci/lactococci, lactobacilli, and enterococci are common microbial constituents of honeybee-collected pollen [28]. Moreover, Lysobacter was recently detected as a minor component of intestine honeybee microbiota [29]. Interestingly, Meiothermus was reported to be a beneficial bacterium of phytophagous insects, i.e., Phasmotaenia lanyuhensis [30].
The major bacterial taxa in fir-oak honey were Lonsdalea, causing acute oak decline [31], and Zymobacter, a facultative anaerobic fermenter [32], which covered 18.53 ± 4.96% and 17.22 ± 4.48% of the total reads, respectively ( Figure 2). Recent findings revealed that "insects visiting drippy blight diseased red oak trees are contaminated with the pathogenic bacterium Lonsdalea quercina" [33]; therefore, this bacterium appears to be transmitted from infected oak trees to honeybees and subsequently to fir/oak-originated honey. The presence of Zymobacter may result in a reduced Lactobacillus population in these honeybees, since this rarely appearing sugar-tolerant bacterium of Halomonadaceae [34] may be favored as a specialized alternative fermenter [35]. Moreover, the acetic acid bacterium Asaia is considered as an emerging symbiont of Apis mellifera [36]. Kocuria spp., which were among the major taxa identified in fir-related honeys (fir and fir-oak honeys; Figure 2), have been identified as the most abundant species in honeybees obtained from beehives in Turkey, which, however, were entomopathogens of honybees [37].
Lysobacter (37.31 ± 5.58% of the total relative abundance), Chryseobacterium (12.66 ± 4.44%), and Paenibacillus (8.81 ± 7.73%) were the predominant microbiota in Arbutus unedo honey ( Figure 2). Although Chryseobacterium is rarely detected in honey [26], this taxon belongs to the native microbiota of Arbutus unedo (strawberry tree), since this bacterium was detected in all strawberry tree specimens examined by Martins et al. [38]. Apart from Lysobacter, which is a representative of honeybee gut microbiota [29] and a plant-associated microbe with plant-protective properties [39], various Paenibacillus spp. have been also detected in the microbiome of honeybee intestine [40]. Although certain members of the genus Paenibacillus are entomopathogenic to honeybees [41], various Bacillus spp. can be beneficial as biological control agents delivered from honeybees to plants [42].
Among other major taxa (Figure 2), Sphingomonas spp., which are also considered as common constituents of honeybee gut microbiota [43,44], were detected in all honey samples examined. Moreover, Methylobacterium has been previously detected in Ceratina bees [45].
Performance of correlation network analysis in honeys examined showed the strong relationships among certain members of Actinobacteria (Kocuria, Nakamurella, Tessaracoccus, Acidothermus, Nocardia, Arcanobacterium, and Propionicicella), Bacteroidota representatives (Flavobacterium and Ferruginibacter), the Rhodobacteraceae genera Amaricoccus and Gemmobacter, and the methylotrophic bacterium Methylotenera. Such microbiota may be involved in the decomposition of complex compounds present in flowering plants and honey, such as phenolics, flavonoids, lignin, and cellulose [46]. A second distinct cluster was formed solely by the major microbiota of the fir-oak honey (Lonsdalea, Zymobacter, and lactic acid bacteria such as Fructobacillus and Leuconostoc), indicating the interaction of Lonsdalea with key fermentative taxa in this honey sample.
Lactobacilli and especially L. kunkeei (basionym of Apilactobacillus kun-keei) strains in the gut of melliferous bees have been reported to induce resistance against deltamethrin, a pesticide that is considered as the major threat for pollinators [47]. Lactobacillus spp. are members of the native gut microbiome of honeybees, which have been found to enhance antiviral properties, even during tetracycline treatment that increases sensitivity to viral infections [48]. Moreover, honey with a high Lactobacillus population was reported to exhibit the greatest antioxidant activity among other honey samples examined by Wu et al. [49]. L. kunkeei (A. kun-keei) in the gut of Apis mellifera has been reported to exert antimicrobial activity and act as probiotic against the honeybee pathogens Paenibacillus larvae and Melissococccus plutonius [50,51]. Indeed, L. kunkeei is favored by a long evolutionary symbiotic relationship in the gut of honeybees [52]. Moreover, Goh et al. [53] reported that Lactobacillus strains, which were isolated from stingless bees and their products, exhibited antimicrobial properties against Listeria monocytogenes. Kaškonienė et al. [54] also reported that lactococci and lactobacilli during solid-state lactic acid fermentation of bee pollen exerted antibacterial activity against Micrococcus luteus, Staphylococcus aureus, and Escherichia coli. Similarly, lactic acid bacteria from the intestine tract of honeybees exhibited inhibitory effects against Bacillus cereus, Escherichia coli, Pseudomonas aeruginosa, Salmonella typhimurium, and Staphylococcus aureus [55].
Comparative analysis of microbial communities in the examined honey samples revealed that their major taxa were either indigenous microbiota of honeybees, mostly beneficial and to a lesser extent potential foodborne and human pathogens, and/or inhabitants of bee pollen and flowering plants, including both plant protective and phytopathogenic microorganisms. In addition, honey originated from cotton (an industrial plant) included a higher proportion of potential foodborne and human pathogens compared to honey originated from bees fed with pollen of wild plants, a fact that may be attributed to the higher abundance of antimicrobial compounds in forest plants [56,57].
A total of 14 taxa were identified in all honey samples (Table 2), independently from their geographic distribution and pollen origin, indicating that these genera may serve as authenticity bioindicators of honey and honeybee-derived products. All of these are either aerobes, capable of growing in low-oxygen conditions (lysobacters and sphingomonads) [58,59] or presenting weak anaerobic growth (Microbacterium), as well as being facultative anaerobes (e.g., bacilli/anoxybacilli, paracocci, pseudomonads, staphylococci, and phenylobacteria) or common fermenters of insect gut (Propionibacterium spp.). Therefore, restricted oxygen levels in the honeybee gut and/or honey seem to play a role in shaping the microbial community structure in honey.  In all honey samples examined, methylotrophs represented an important fraction of microbial population. In particular, methylotrophs covered 1.3-3.0% of the total relative abundance ( Table 3), independently of the geographical and pollen origin. Methanol has been reported to be among the most abundant volatile compounds of honey [60,61]. Interestingly, plant-colonizing methylotrophs of the genus Methylobacterium were capable of delivering insecticidal proteins and have been associated with the plant protective properties of the common biocontrol agent Bacillus thuringiensis [62]. Moreover, volatiles such as methanol promote the growth of the phyllosphere microbiota [63], serving as an energy source for methylotrophic epiphytes [64]. Thus, both methanol concentration and co-evolution of honeybees with their plant hosts and phyllosphere microbiota appear to shape methylotrophic communities in honey. In conclusion, native microbiota of the honeybee microbiome, such as fermentative bacteria, either with probiotic properties (e.g., lactobacilli) or potential foodborne and human pathogens (e.g., staphylococci and streptococci), as well as common inhabitants of bee pollen (e.g., pollen wall decomposers, such as Pseudomonas spp.) and flowering plants (e.g., beneficial microorganisms of plants and phytophagous insects, such as Lysobacter and Meiothermus spp., respectively), including plant pathogens of honeybee hosts that are transmitted from flowering plants to honey via honeybees (e.g., Lonsdalea, the causing agent of acute oak decline), constitute the microbial community structure in natural honeys.

Collection of Honey Samples
Honey samples from various locations in Greece were obtained aseptically and further examined in terms of their antimicrobial properties. Four honey samples from beehives situated in different geographical areas of Greece, i.e., from the Prefecture of Epirus (the one honey sample was produced from honeybees fed with fir and oak and the other with Arbutus unedo pollen and nectar), the Regional unit of Phthiotis, Mendenitsa region (from fir pollen and nectar) and the Regional unit of Karditsa, Palama region (from cotton pollen and nectar), were collected and further studied, due to their high antimicrobial activity (as examined in Stavropoulou et al. [65,66]). These honey samples were subsequently subjected to DNA extraction and phylogenetic analysis of microbial communities through Illumina sequencing.

DNA Extraction from Honey Samples of Different Pollen Origin and Performance of Amplicon Sequencing
The collected honey samples exhibiting high antimicrobial activity against foodborne and human pathogens were subjected to DNA extraction. Genomic DNA was extracted from honey through the use of NucleoSpin Tissue Kit (Macherey-Nagel, Düren, Germany), by following the instructions of the manufacturer. Aliquots of 60 µL of 10 mg/mL lysozyme and 60 µL of 10 mg/mL lysostaphin, as well as 6 µL of 60 U/µL lyticase (enzymes supplied by Sigma-Aldrich, Germany), were added per treated sample to facilitate the lysis of Gram(+) bacteria and yeast strains, respectively. The V4-V5 region of the 16S rRNA gene was amplified using the primers 515F (5 -GTG YCA GCM GCC GCG GTA A-3 ) and 909R (5 -CCC CGY CAA TTC MTT TRA GT-3 ). No amplification was achieved using primer set ITS1F (5 -CTT GGT CAT TTA GAG GAA GTA A-3 ) and ITS4R (5 -TCC TCC  GCT TAT TGA TAT GC-3 ) for fungi, including yeasts. Amplification of partial 16S rRNA gene was conducted through a thermal scheme comprising of 3 min denaturation at 94 • C, succeeded by 28 cycles of 30 s denaturation at 94 • C, 40 s primer annealing at 53 • C, and 1 min DNA elongation at 72 • C, before a 5 min DNA thermal extension at 72 • C to complete the thermocycling reaction. The amplification reaction for the ITS region was performed by using 2 min denaturation at 94 • C, succeeded by 35 cycles of 30 s denaturation at 95 • C, 30 s primer annealing at 55 • C, and 1 min DNA elongation at 72 • C, before a 10 min DNA thermal extension at 72 • C to complete the thermocycling reaction. Illumina sequencing reactions were carried out at Mr DNA' (USA) in MiSeq equipment (Illumina, Inc., San Diego, CA, USA) after amplicon cleanup through the use of DNA purification beads.

Bioinformatic Analysis and Analysis of Variance (ANOVA)
Bacterial amplicons were proceeded to demultiplexing and trimming, and partial 16S rRNA gene sequences of abnormal size and low-quality were removed. The assembled reads were further subjected to denoising and chimera discard using USEARCH v.11 [67,68]. Bacterial sequences were clustered and ZOTUs (zero-radius OTU of denoised sequences) were generated, following National Center for Biotechnology Information (NCBI) taxonomy. Ecological indices, i.e., Chao1, Shannon, Simpson, and Fisher diversity indices, were sequentially calculated using the MicrobiomeAnalyst online suite of omics tools [69]. Relationships among bacterial taxa were identified in the MicrobiomeAnalyst platform via multiple correlation network matrix analyses based of the SparCC score at 0.97 correlation coefficient. The sequenced amplicons were deposited in the Sequence Read Archive (SRA) of the NCBI platform under the BioProject accession number PRJNA913297.
Analysis of variance (ANOVA) was conducted using Past v.4.10 [70] to identify statistically significant differences among the relative abundances and diversity indices of the bacterial communities identified in these Greek honey samples exhibiting high antimicrobial activity against foodborne and human pathogens.