Gut Microbiome of Two Different Honeybee Workers Subspecies in Saudi Arabia

1Department of Biological Sciences, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia.71491. 2Department of Biology, College of Science, Tabuk University, Tabuk, Saudi Arabia.74191. 3Department of Genetics, Faculty of Agriculture, Ain Shams University, Cairo, Egypt. 11241. 4Princess Al JawharaAlbrahim Centre of Excellence in Research of Hereditary Disorders (PACER-HD), King Abdulaziz University, Jeddah, Saudi Arabia.

Honeybees belong to the genus Apis, which is known for its tremendous role in pollination. Unfortunately, honeybee population is recently declining with a potential risk on the agricultural service and subsequently the food supply, not only locally in Saudi but also globally 1 . There is a known mutually beneficial relationship between honeybee gut microbiome and its host. The host provides the optimum environment for bacterial growth, while the bacterial community in honeybee guts aids in efficacy of nutrients absorption, optimum growth and development of the host and its ability to defend pathogens, and its adaptation to surrounding environment 2 .Honeybee gut represents a simple model system to study the relationship between gut microbiome with honeybeehosts 3,4 .The bacterial community in adult honeybee workers is diverse and estimated to reach one billion bacterial cells in each worker's gut 5,6 . Such a diversity in bacterial community is dependent on the type of flower that hosts the insect, as well as many other environmental factors 7 .Gut microbiome of honeybee (Apis mellifera) workers is composed of eight to nine core species 8,9 , e.g., Bartonella apis 10 , Acetobacteraceae 11 ,Parasaccharibacter 11 , Snodgrassella alvi 12 , Bifidobacterium asteroids 13 , Lactobacillus sp. 14 , Frischella perrara 15 and Gilliamella apicola 12 .
The two most common bee species that are widely distributed throughout the kingdom of Saudi Arabia are the indigenous Apis mellifera jementica, which is a native species, and Apis mellifera carnica, which is imported from Egypt 16 as honey production of domestic bees does not meet the growing demands in Saudi Arabia. Moreover, the production cost is relatively high. Exotic bee colonies have been imported over time, reaching 200,000 bee packages annually 16 . It is well known for local beekeepers that the indigenous bees A.m. jementica highly tolerates local stressful conditions when compared with exogenous races A.m. carnica, particularly during summer when the air temperature becomes extremely high. It is also noticed that at high temperatures, indigenous bees continue to forage for pollen and collect nectar, whereas imported bees will stop foraging 16 . Initial reports revealed that the subspecies of exotic honeybees have lower heat tolerance, shorter foraging durations and are more susceptible to Varroa mites when compared with indigenous bees 16 .
In the present study, we compered the gut microbiome composition and diversity of the adult honeybees of Apis mellifera jementica and Apis mellifera carnica in Saudi Arabia using highthroughput 16S rRNA gene sequencing technology.

Material and Methods sample collection,isolation of guts microbiota and dna extraction
Five samples each from honeybee workers of A.m. jemenitica and A.m. carnica were collected in November 2019 from a single hive of Beekeeper Cooperative Association at Al Baha, Saudi Arabia. The collected samples were immediately stored at "80°C. For whole gut dissection of honeybee workers,surface disinfection was done using 1 ml aqueous ethanol (70%, v/v) for 45 sec. Dissected guts were, then,placed in a pre-frozen mortar and 700ìl S1 lysis buffer (Invitrogen, Thermo Fischer Scientific, USA) were added and guts were transferred to bead tube for extraction process. DNA of gut samples was extracted by the genomic DNA extraction kit (Invitrogen, Thermo Fischer Scientific, USA), and stored at -20°C for further molecular analysis.

PCR amplification
PCR was run to amplify bacterial 16S rRNA gene of the variable regions V3-V4. The two universal primers used for PCR are 341F 52 -ACTCCTACGGGAGGCAGCAG-32 ( f o r w a r d p r i m e r ) a n d 8 0 6 R 5 2 -GGACTACHVGGGTWTCTAAT-32 (reveres primer). The PCR conditions were set as the following: one cycle for initial denaturation at 95°C for 5 min; 25 cycles of denaturation at 95°C for 30sec followed by annealing at 56°C for 30 sec and primer extension at 72°C for 40 sec; and a one cycle for final extension at 72°C for 10 min. The generated PCR products were checked for quality and selected products were utilized in preparing Illumina DNA libraries. DNA sequencing was run using Illumina Miseq platform (Illumina, San Diego, CA) at Beijing Genome Institute (BGI), China to generate high-quality pair-ends of ~300 bp.

statistical analysis
The high quality paired reads produced in fasta files as raw data were de-multiplexed, qualityfiltered and trimmed by trimmomatic package (Version 0.33) through Quantitative Insights Into Microbial Ecology 2 pipeline (QIIME2, v1.80). Obtained reads were merged into single sequence files by the Fast Length Adjustment of SHort reads (FLASH, Version 1.2.11). In order to assign generated unique sequences into operational taxonomic units (OTUs), reads were tagged and clustered into OTUs with similarity cut off of 97% using the de novo OTU piking procedure. Usearch (Version 7.0.1090) 19 was, then, used to remove Chimeric sequences. Taxonomies were plotted against the gut Microbiome Database (HOMD RefSeq, Version 13.2) through the RDP classifier (Version 2.2) 17 and the Green-genes database (version 201305 18 16S rDNA database, http://qiime. org/home_static/dataFiles.html) with a cut off of 70%. Alpha diversity indeces were measured in order to assess the intra-species variations within a given sample using Mothur (v1.31.2). Alpha diversity and rarefaction curve boxplots were constructed using software R (v3.1.1). To investigate the inter-species variations within samples, the beta diversity matrices were conducted and visualized using principal coordinate analysis (PCoA) by package 'ade4' of software R (v3.1.1). Also, heat maps were generated using the package 'gplots' of software R (v3.1.1), and, then, sequence alignments were searched against the Silva core set (Silva_108_core_aligned_seqs) by using PyNAST 'align_seqs.py'. The obtained OTU phylogenetic tree was, then, plotted by software R (v3.1.1), and visualized through QIIME2 (v1.80).
Annotation of generated OTUs was done in order to detect the relative abundance at different taxonomical levels (phylum, genus and species). Finally, Metastats, PERMANOVA and Benjamini-Hochberg false discovery rate (FDR) correction were also used to correct for multiple hypothesis. The Linear Discriminant Analysis (LDA) Effect Size (LEfSe) was applied using software LEfSe with the online interface Galaxy (version 1.0.0; http://huttenhower.sph.harvard.edu/galaxy/root),to discriminate the two taxonomic races determining highly presented bacterial taxon within each race depending on statistical significance.

statistics of 16s rrna sequence data
The five gut microbiome samples of A.m. carnica were identified asC1 to C5, while the five gut microbiome samples of A.m. jemenitica were identified as J1 to J5. Illumina MiSeq was used in sequencing thepartial 16S rRNA gene   Figure S1). The tagged sequences were assigned to a total of 171 OTUsacross samples ranging from 45 (J3) to 154 (C5) OTUs (Table S1) Table 2).
Principal coordinate analysis (PCoA) was used to display the diversity as well as the differences in OTU composition.Diversity of A.M.C subjects was higher towards positive and negative PCA 1 directions (PC1), where as that of A.M.J subjects was higher towards positive and negative PCA 2 directions (PC2). As an overall picture, the diagram shows that the mean value of A.M.C group was localized in positive portion of PC1 and negative portion of PC2, whereas A.M.J group was mainly localized in the positive portion of PC2( Figure 2). The principal coordinate analysis (PCoA) plots were created using a Bray-Curtis distance matrix and the samples were plotted to represent the microbial community compositional differences between samples. The plots are dimensionally scattered in accordance to their gut microbiome compositional relationships.
The results of the present study indicate that the differences ingut microbiomes between these two groups are possibly due to the different origins of worker honeybees of the two subspecies. The stacked number of OTUs and the number of observed species for different samples as rarefaction measures are shown in Figure S2. When the refraction curves inclines ( Figure S2a) or stops climbing ( Figure S2b), the produced data would be enough for further analysis. However, as long as the curve is still climbing, the complexity of the data in samples become higher; since more species being detected throughout sequencing analysis. The two rarefaction curve measures refer to the maximum number of sequences attained for all samples that allows to study taxonomic relative abundance and to assess eligibility of such data to represent all species of any microbial community. The findings from both rarefaction measures show that 54,000 is the maximum number of sequence reads that can be used further in studying taxonomic abundance ( Figure S2).

structure of gut microbiomes across the two honeybee workers
Two taxonomic ranks (phylum and species) were used in the comparison of gut microbiomes between adult honeybee workers A.M.C and A.M.J at the phylogenetic level ( Figure 3). The results indicate that phylum Firmicutes harbours 24 genera,while Proteobacteria,Actinobacteria, Bacteroidetes and Thermiharbour 23, 8, 6 and 2

differential abundance of microbes due to different origin of worker
The observed microbial taxa along with their redundancies across different samples identified after OTU annotation are described in Table S2. The taxa refer to phylum, class, order, family, genus, and species. Eight phyla of the gut bacteria were identified according to relative abundance. They are Actinobacteria, Bacteroidetes, Cyanobacteria, Firmicutes, Protobacteria, TM7, Tenericutes and Thermi (Figure 4). Aligning with the number of genera of each phylum shown in Figure 3, the most abundant phylum were Firmicutes (57%), Protobacteria (31%) and Actinobacteria (10%) inA.M.C group (Figure 4). Meanwhile, Firmicutes (48%), Protobacteria (44%) and Actinobacteria (6%) were the most abundant in A.M.J group (Figure 4). The comparison at phylum level revealed a significant increase in Cyanobacteria in the A.M.C group (P-value = 0.031746), while a significant increase of Protobacteria in the A.M.J group (P-value = 0.037724) (Table S3). Interestingly, Table S3 also indicates the existence of the three phyla TM7, Tenericutes and Thermi only in A.M.C group. The previous results align with those of the heat map at phylum level as Firmicutes, Protobacteria and Actinobacteria were shown to be the most abundant phyla across samples and groups ( Figure S3). In terms of species relative abundance in the gut microbiomes of two groups A.M.C and A.M. J Bacteroides_fragilis, Bacteroides_  o v a t u s , C o m m e n s a l i b a c t e r _ i n t e s t i n i , Blautia_producta, Melissococcus_plutonius, Ruminococcus_gnavus,Saccharibacter_floricola and Snodgrassella_alviwere shown to be the most abundant ( Figure 5).The figure also indicates that a large proportion of the OTUs were not assigned to a certain species (93.80% for A.M.C and 86.20% for A.M.J). We have no explanation for these results except that a large number of species in workers of honeybee was not identified or classified before. The results in Table S4 indicates a significant increase of Melissococcus_ plutoniusin the gut microbiome of A.M.C (P-value = 0.034454), while Snodgrassella_alvi in theA.M.J group (P-value =0.008948). Results for the latter species Snodgrassella_alvi align with that presented in Figure 5c. TheRuminococcus_ gnavusandSaccharibacter_floricolawere not existed in theA.M.J group. The heat map at species level indicates that Snodgrassella_alviharbours the highest relative abundancea cross all samples ( Figures S4).
Linear discriminant analysis effect size (LEfSe) and its LDA scores (Ã 3) were used to identify possible biomarkers in gut microbiota that refer to the origin of the host (  while Betaproteobacteria (Neisseriales, Neisseriaceae,Snodgrassella sp. and Snodgrassella_ alvi)in A.M.J (Figure 6b).

discussion
The gut microbiome structure of honeybee workersis dependent upon monophyletic origin of  This genus produces several compounds in honeybee gut with known antimicrobial activities such as organic acids, hydrogen peroxide, bacteriocin, reutericyclin and reuterin that mostly inhibit decaying and protects against pathogenic bacteria, as well as some fungi 25,26 .Therefore,Honeybees likely use lactobacilli as probiotic 27 .In the present study, the dominance of Lactobacilli in both A.m. carnica and A.m. jemenitica adult workers is supported by the presence of low pH (3.9) of honey and nectar 28 . This is concluded because of the ability of lactobacilli to ferment sugar in the gut of honeybee workers and, hence,to generate acidic environment 29 , which inhibits the growth of many other bacteria. The low abundance in Lactobacillaceae was reported to be associated with the presence of pathogenic bacteria 30 .
Genus Bifidobacterium,gram-positive bacteria belonging to the Actinobacteria phylum, was also identified in gut of both A.m. carnica and A.m. jemenitica adult workers. Again, it is dominant in rectum, and a core gut bacteria of honeybee workers. Bifidobacterium strains carry large surface proteins, which have a role in adhesion or degradation of plant materials 7,31,32 . Additionally, Bifidobacterium carriesgene clusters that are responsible for the production and utilization of trehalose, which is a disaccharide molecule used by insects as an energy reservoir, in comparison to glycogen, which is the energy storage form in mammals 33 .
Family Neisseriaceae and its descendent Snodgrassella_alvi(S. alvi), gram-negative bacteria belonging to Betaproteo bacteria phylum, significantly increased in A.m. jemenitica. These bacteria participate in oxidation of carbohydrates. However, the pathway for the uptake and glycolytic breakdown of carbohydrates does not exist in S. alvi, thus,this bacteriumis located consistently within the periphery of the insect's gut lumen. This area has high oxygen concentrations and this environment is preferable for S. alvidue to its dependence on aerobic respiration 34, 35 . Insects depend on the aerobic oxidation of carboxylates rather than breaking down carbohydrates resulting in various products such as citrate, malate, acetate and lactic acid that serve as energy sources 12,27 . The steady co-exits of S. alvi with other fermentative bacterial taxa in the same gastrointestinal environment can result from utilizing separate sets of resources leading to metabolic variations suggesting a syntrophic interaction. For example, S. alvican utilize some of the substrates such as lactic acid, acetate and formate, which are produced from carbohydrate fermentation 36,37 . Furthermore, S. alvi and G. apicola 38 are enriched with genes encoding biofilm formation. The two species inhabit the host's ileum, indicating that the biofilm can provide a protective layer against pathogens.
The bacteria of the family Acetobacteraceae and its descendent genus Commensalibacter(also referred to as Alpha 2.1), gram-negative bacteria belonging to phylum Proteobacteria,were identified as a core member of the gut microbiota in honeybees and bumble bees 9,31 . It was observed mainly in the midgut and hindgut of honeybee workers. In our study, Commensalibacter presents in A.m. carnica and A.m. jemenitica. However, Saccharibacterflorica (Alpha-2.2) presents only in A.m. carnica. Furthermore, Saccharibacterflorica is isolated from pollen, suggesting that this phylotype is associated with flowers 39 .The role of these phylotypes (Alpha 2.1 andAlpha-2.2) is associated with their abilities to adapt with fast growing metabolic processes, with two distinctive mechanisms. Alpha2.1 bacteria harvest energy through a wide range of substrates linked and utilized through a flexible oxidative and biosynthetic metabolism. Whereas, Alpha2.2 bacteria, that lack alternative oxidative pathways, determine metabolic processes through oxidative fermentation after harvesting glucose for rapid energy 40 .
T h e b a c t e r i a o f t h e f a m i l y Enterococccaceae and its descendent species Melissococcusplutonius, gram-positive bacteria of phylum Firmicutes, present in low abundance (3%) in gut microbiome of A.m. carnica honeybee workers. This conclusion was also noted in previous reports 41 .M. plutoniusis known to cause the European foulbrood (EFB) in earlystage of honeybee larvae, with assistance from secondary invaders (Enterococcus faecalis, Paenibacillus alvei and Bacillus pumilus). M. plutonius was shown to have 30 different sequence types clustered under three clonal complexes (CC 3, CC12, and CC13) 42,44 , where CC13is the least virulent complex 43,45 . Honeybee workers transmit M. plutonius between colonies via robbing and drifting 46,47 .Erban et al. 45 compared control samples from the EFB zone with samples from EFB zone without clinical symptoms,and bees from colonies from EFB zone with clinical symptoms. The study identified a 100-fold higher prevalence of M. plutonius in colonies with EFB symptoms, while it only presents in 3 of 16 control colonies that are distant from the EFB zone. This suggests that M. plutonius has lower abundance in healthy honeybee colonies, which is consistent with the results of the present study.

conclusion
The present findings indicative that differences in gut microbiome structures of honeybee workers of the two subspecies A.m. carnica and A.m. jemenitica are due to varied monophyletic origin of the host. These findings support previous results suggesting that honeybee workers have a mutual coevolving relationship with specific group of bacteria. This group of bacteria co-exists and is maintained throughout the descending generations of the host. Inclusion of more subspecies inhabited in Saudi Arabia along with ones of this study can further support our findings.

acknowledgeMents
This study was supported by Beekeeper