Colonization of the gut microbiota of honey bee (Apis mellifera) workers at different developmental stages

Journal Pre-proof Colonization of the Gut Microbiota of Honey Bee (Apis mellifera) Workers at Different Developmental Stages Zhi-Xiang Dong, Huan-Yuan Li, Yi-Fei Chen, Feng Wang, Xian-Yu Deng, Lian-Bing Lin, Qi-Lin Zhang, Ji-Lian Li, Jun Guo

PII:

S0944-5013(19)30715-3

DOI:

https://doi.org/10.1016/j.micres.2019.126370

Reference:

MICRES 126370

To appear in:

Microbiological Research

Received Date:

4 July 2019

Revised Date:

28 September 2019

Accepted Date:

26 October 2019

Please cite this article as: Dong Z-Xiang, Li H-Yuan, Chen Y-Fei, Wang F, Deng X-Yu, Lin L-Bing, Zhang Q-Lin, Li J-Lian, Guo J, Colonization of the Gut Microbiota of Honey Bee (Apis mellifera) Workers at Different Developmental Stages, Microbiological Research (2019), doi: https://doi.org/10.1016/j.micres.2019.126370

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Colonization of the Gut Microbiota of Honey Bee (Apis mellifera) Workers at Different Developmental Stages

Zhi-Xiang DONG1, Huan-Yuan LI1, Yi-Fei CHEN1, Feng-WANG, Xian-Yu DENG, Lian-Bing LIN, Qi-Lin ZHANG1* [email protected], Ji-Lian LI2* [email protected],

-p

Faculty of Life Science and Technology, Kunming University of Science and Technology,

Kunming, Yunnan 650500, China

Institute of Apiculture, Chinese Academy of Agricultural Science/Key Laboratory of

lP

2

re

1

ro of

Jun GUO1,* [email protected]

Corresponding authors: (J. G.), (J. L. L.) and (Q. L. Z.)

ur

*

na

Pollinating Insect Biology, Ministry of Agriculture, Beijing 100093, China

Jo

Abstract：The role of the gut microbiome in animal health has become increasingly evident. Although the structure of the gut microbiome of A. mellifera is well known, little is known about the dynamic change across different developmental stages. In this study, we explored the dynamic changes of the gut microbiota of A. mellifera at different developmental stages covering the whole life cycle using high-throughput 16S rRNA gene sequencing. The results indicated that the core (shared) gut microbiota changes significantly among different developmental stages. The diversity

of the bacterial community in workers among different ages was significantly different. In addition, by comparing the core gut microbiota among different-aged workers, we found that newly emerged workers had fewer core microbiota. Three genera, Gilliamella, Frischella, and Snodgrassella, were significantly colonized at 1 day poste mergence (dpe); Lactobacillus, Bifidobacterium, Commensalibacter were significantly colonized at 3 dpe and significantly reduced with Gilliamella. Lactobacillus kunkeei and Bartonella were significantly colonized at 12

ro of

dpe and were significantly decreased with Lactobacillus helsingborgensis. Commensalibacter and

Bifidobacterium were significantly decreased at 25 dpe, and Bacteroides, Escherichia-Shigella, and Porphyromonadaceae were significantly decreased between 19 and 25 dpe. Our results reveal the

-p

succession of the gut microbiota of workers from birth to senescence, which provides a theoretical

lP

re

basis for further exploring the roles of gut microbiota during different developmental stages.

Keywords: Apis mellifera workers, Colonization, Bifidobacterium, Lactobacillus kunkeei, Core

ur

na

microbiota

1 Introduction

Jo

The honey bee, A. mellifera, is an important economic insect. The U.S. government estimates the annual social gains of honey bees to be between 1.6 billion and 5.7 billion dollars (Southwick et al., 1992). Honey bees are also important pollinators; from 2006 to 2008, the average annual value of bee pollination reached 304.22 billion RMB, which is equivalent to 12% of the total agricultural output value of China (Liu PD, 2011). However, in recent years, the survival of honey

bees has suffered a great threat. Particularly, the US reported Colony Collapse Disorder (CCD), in which thousands of hives appeared empty with a mysterious disappearance of bees (Stokstad, 2007). This accident is caused by many factors, including climate change (Langowska et al., 2017), the heavy use of pesticides (Colin et al., 2019), the abuse of antibiotics (Owen, 2017), and pathogen infection (Cameron et al., 2011). All these factors affect the feeding and other behaviors of honey bees. Honey bee gut microbiota have a close relationship with the host and provide

ro of

several advantages to the host, such as promoting the digestion of food, essential nutrients, the degradation of toxic components, pathogen defense, and the regulation of host development, behavior, and immunity (Moran NA, 2003).

-p

According to the 16S rDNA analysis of the community in the gut, researchers found that

re

there were mainly nine types of bacteria distributed in the intestine of workers (Martinson et al.,

lP

2011; Sabree et al., 2012), five of which were present in almost all honey bees, including two gram-negative species, Snodgrassella alvi and Gilliamella apicola (Kwong et al., 2013). Among

na

these bacteria, Lactobacillus firm-4 and firm-5, belonging to the phylum firmicutes, were the most dominant and widespread in the gut (Babendreier D, 2007). The abundance of Bifidobacterium is

ur

relatively lower in comparison with the majority of bacterial species, but it is widespread (Bottacini et al., 2012a). Compared with the core bacteria, the abundance and distribution of

Jo

Frischella perrara (Philipp et al., 2013), Bartonella apis (Jeyaprakash et al., 2003), Parasaccharibacter apium (Vanessa et al., 2014), and "Alpha 2.1" are not stable. Honey bees are a good model for investigating the gut microbiota. Anderson et al. (2018) studied workers and queens with different age phenotypes and found that despite workers and queens sharing many gut bacterial species, the structure of gut microbiota is markedly different between them (Anderson et

al., 2018). Their research also revealed that carbon accumulation in queens is significantly correlated with the increase in Lactobacillus and Bifidobacterium and the decrease in the phylum Proteobacteria (Anderson et al., 2018). In addition, it has been indirectly proven that the gut flora of bees also play an important role against pathogenic microorganisms. Previous studies showed that the abundance of these core gut bacteria were relatively lower than Crithidia bombi in bumble bees, and fecal transplants of newly developed adult bumblebees were found to have a similar

ro of

effect (Hauke et al., 2011). However, when the structure of the gut microflora is disturbed by

antibiotics, the diversity of the gut flora decreases and honey bees with abnormal intestine bacteria

will be more susceptible to opportunistic pathogens, reducing the survival rate (Raymann K, 2017).

-p

Similar studies have been conducted on Apis cerana by studying A. cerana infected with Nosema

re

ceranae. The structure and diversity of gut flora were significantly changed, with an increase in

lP

mortality (Huang SK, 2018). These results demonstrate the importance of honey bee gut bacteria to host health.

na

The development of high-throughput sequencing technology has greatly promoted research on gut microbiota (Ross, 2012). A clearer understanding of the composition of bee gut flora has

ur

been gradually formed using independent-cultured methods, particularly high-throughput sequencing of the 16S rRNA gene (Kwong et al., 2016). For example, Jones et al. sequenced the

Jo

16S rRNA gene of workers exposed to different environmental factors and found that exposure to different environmental factors influenced the relative abundance of partially microbial groups (Jones et al., 2017). Kakumanu et al. found that pesticide exposure significantly affected the structure of gut flora in honey bee workers using high-throughput sequencing technology (Kakumanu et al., 2016). Next, Motta et al. explored the effect of glyphosate on bee gut flora and

found that exposure can perturb beneficial gut bacteria (Motta et al., 2018). Despite many previous studies focusing on the gut microbiota of A. mellifera, the dynamic succession of the gut flora of A. mellifera across the whole life cycle is poorly understood. Here, high-throughput sequencing technology was used to analyze the gut microbiota of workers at 0 day postemergence (dpe), 1 dpe, 3 dpe, 7 dpe, 12 dpe, 19 dpe, 25 dpe, 30 dpe, 35 dpe and 40 dpe. The dynamic changes in the gut microbiota of workers were illuminated, aiming to explore the

ro of

succession of workers from birth to senescence and to further understand the dynamic changes in

the community of the gut microbiota of workers. These results provide a valuable genetic resource

-p

for the use of gut microbiota to improve the health conditions of bees in the future.

re

2 Materials and Methods

lP

2.1 Worker sampling

We collected worker bees (Apis mellifera) in Kunming city, Yunnan Province, China, from early

na

July to late August 2018. To confirm the age of the workers, we chose a frame from a group of healthy colonies, removed the workers from the frame, and placed the workers at 34 °C in the dark

ur

(mimicking hive conditions). First, we randomly took out 6 workers that were just out of the cocoons as 0 dpe, and then immediately placed them into a 1.5 mL centrifuge tube and stored

Jo

them in an ultralow temperature freezer (EU1DW/BD-55W321EU1, China) at -80 °C (samples did not come into contact with the nurse bees). After that, the frame was returned to the dark incubation environment. Next, 300 newly emerged workers were marked with Testors enamel paint and then put back into the original hive for natural growth. Six workers were randomly collected at 1 dpe, 3 dpe, 7 dpe, 12 dpe, 19 dpe, 25 dpe, 30 dpe, 35 dpe and 40 dpe. After

dissection of the gut in a sterile environment, the entire gut was removed and placed into a 1.5 ml centrifuge tube. All the samples were stored in 75% alcohol in the ultralow temperature freezer at - 80 °C for DNA extraction.

2.2 DNA extraction and PCR amplification Next, the DNA of the gut microbiota was extracted using an E.Z.N.A.® soil DNA Kit (Omega,

ro of

U.S.) according to the manufacturer’s protocols. The final DNA concentration and purity were determined with a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, USA), and

DNA quality was checked by 1% agarose gel electrophoresis. The V3-V4 hypervariable regions of

-p

the bacterial 16S rRNA gene were amplified with primers 338F (5’-

re

ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’) with a

lP

thermocycler PCR system (Applied Biosystems, USA). The PCR were implemented following the program: 3 min of denaturation at 95 °C, 27 cycles of 30 s at 95 °C, 30 s for annealing at 55 °C,

na

and 45 s for elongation at 72 °C, and a final extension at 72 °C for 10 min. PCR were performed in triplicate in 20 μL mixtures containing 4 μL of 5 × FastPfu Buffer, 2 μL of 2.5 mM dNTPs, 0.8

ur

μL of each primer (5 μM), 0.4 μL of FastPfu Polymerase and 10 ng of template DNA. The PCR products were purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, USA)

Jo

and then quantified using QuantiFluor™-ST kits (Promega, USA) according to the manufacturer’s protocol.

2.3 Illumina MiSeq sequencing and processing of sequencing data Purified amplicons were pooled in equimolar and paired-end sequenced (2 × 300 bp) on an

Illumina MiSeq platform (Illumina, USA) according to the standard protocols from Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). Raw fastq files were quality-filtered by Trimmomatic (version 0.36) and merged by FLASH (version 1.2.7) (Zhang et al., 2018) with the following criteria: (i) The reads were truncated at any site receiving an average quality score <20 over a 50 bp sliding window. (ii) Sequences with overlaps longer than 10 bp were merged according to their overlap with mismatch no more than 2

ro of

bp. (iii) Sequences of each sample were separated according to barcodes (exactly matching) and primers (allowing 2 nucleotide mismatching), and reads containing ambiguous bases were

removed. Then, operational taxonomic unit (OTU) clustering was performed on the sequence with

-p

the following criteria: (i) Nonrepeating sequences were extracted from the optimized sequences to

re

reduce the redundant computation in the process of analysis

lP

(http://drive5.com/usearch/manual/dereplication.html).(ii) Single sequences that did not repeat were removed (http://drive5.com/usearch/manual/singletons.html). (iii) Operational taxonomic

na

units (OTUs) were clustered with 97% similarity cutoff using UPARSE (version 7.1 http://drive5.com/uparse/) with a novel ‘greedy’ algorithm that performed chimera filtering and

ur

OTU clustering simultaneously. The taxonomy of each 16S rRNA gene sequence was analyzed by the RDP Classifier algorithm (http://rdp.cme.msu.edu/) against the SILVA (SSU123) 16S rRNA

Jo

database using a confidence threshold of 70%. In addition, rarefaction analysis was carried out by Mothur.

2.4 Statistical analyses and comparison of microbial communities The Kruskal-Wallis test was used to detect whether there were significant differences in the alpha

diversity metric between the ten groups (ACE, Chao1, Shannon and Simpson). According to the community abundance data obtained, the Kruskal-Wallis statistical method was used to detect the relative richness differences of OTUs in the microbial communities of different groups, and the hypothetical test can be conducted to evaluate the significance (P-values). In addition, FDR correction was performed for the P value of multiple tests, and the statistical Mann-Whitney U test

< 0.05 between groups were declared statistically significant.

3 Results

-p

3.1 Summary of sequencing data

ro of

was used to compare the two groups. The software was the stats package in R. The results with P

re

A total of 3,204,891 high-quality sequences were generated from 60 worker samples, and the

lP

average sequence length was 444 bp. Among the sequences, 3,217 bacterial operational taxonomic units (OTUs) were identified at the 97% sequence similarity cut-off. We estimated the OTU% in

na

60 honey bee samples by the Good's coverage index, and the average bacterial coverage was 0.99 ± 0.00069. This indicated that the obtained data could adequately cover the gut flora of honey bees.

ur

In addition, the rarefaction curves tended to approach the saturation plateau, which indicated that the sequencing data volume was reasonable; more data would only produce a small number of

Jo

new OTUs. Not only that, rarefaction curves also showed that the community richness at 0 dpe was the highest among all sampled time points (Supplementary Figure S1).

3.2 Comparison of bacterial community diversity indices The richness of the community was evaluated by the ACE and Chao1 indices, and the community

diversity was determined by the Shannon and Simpson indices for ten groups. By the Kruskal-Wallis test for comparison, we found that the Shannon, ACE and Chao indices were statistically significant different among different-aged workers (Table 1). This result indicated that the diversity and abundance of gut flora were affected by different ages of workers (Figure 1). The P-value and FDR of diversity index significance were as follows: ACE, P = 0.023, FDR = 0.023; Chao, P = 0.0024, FDR = 0.0072; Shannon, P = 0.0016, FDR = 0.0063; and Simpson, P = 0.010,

ro of

FDR = 0.020. In addition, significant differences in microbial OTU number among different-aged workers were detected (Number of OTUs, P = 0.000085, FDR = 0.00043), indicating that the

-p

change in OTU number was significantly different among host-age stages (Table 1).

re

3.3 Dynamics in microbial flora among different-aged workers

lP

All sequences obtained in this study were classified at the phylum and genus levels. A total of 41 phyla and 838 genera were identified (Figure 2 and Figure 3). The representative sequences at the

na

phylum level are presented in Figure 2. The phyla Proteobacteria, Firmicutes, Actinobacteria, Bacteroidetes, Cyanobacteria, Spirochaetae, Fusobacteria, Deinococcus-Thermus,

ur

Verrucomicrobia, Chloroflexi, Acidobacteria, and Planctomycetes were the dominant phyla, which is similar to previous research (Yun et al., 2018). The relative abundance of phyla in the gut

Jo

flora was significantly influenced by different-aged stages of workers. For the newly emerged workers, Proteobacteria, Firmicutes, Bacteroidetes, Cyanobacteria and Acidobacteria were the dominant phyla, and the relative abundance of Proteobacteria and Firmicutes increased significantly with worker development (Figure 4). The P-value and FDR values were determined by an intergroup comparison and were as follows: Proteobacteria (p = 0.0000341, FDR =

0.000129), Firmicutes (p = 0.0000578, FDR = 0.000171), Actinobacteria (p = 0.0000345, FDR = 0.000129), Bacteroidetes (p = 0.00315, FDR = 0.00561), Cyanobacteria (p = 0.00000321, FDR = 0.0000284), Spirochaetae (p = 0.000314, FDR = 0.000715), Fusobacteria (p = 0.00189, FDR = 0.00352) and Acidobacteria (p = 0.000214, FDR = 0.000517). The 0 dpe bees exhibited the highest bacterial diversity, among which Bacteroidetes (15.82%), Cyanobacteria (7.79%), and Verrucomicrobia (1.94%) were the dominant phyla

ro of

compared with the other day-old workers (Supplementary Figure S2). At 1 dpe, Bacteroidetes

decreased to 12.37%; at 3 dpe, the Bacteroidetes tended to stabilize. Before 19 dpe, all groups showed a similar distribution of the gut flora composition at the phylum level, but the relative

-p

abundance of dominant bacteria varied with worker growth. At 25 dpe, an increase in gut flora

re

diversity at the phylum level was detected, but the composition was highly similar to that before

lP

19 dpe.

At 0 dpe, the gut microbiota of workers mainly consisted of Acinetobacter, Bacteroides and

na

Cyanobacteria, and Acinetobacter was replaced with other bacteria at 1 dpe. These three phyla resided in the gut microbiota of workers before 19 dpe, and their relative abundance was then less

ur

than 0.9% at 25 dpe. In addition, it is worth noting that at 0 dpe, the composition of gut microflora of workers lacked the five core bacteria (A. mellifera). We found that the abundances of the core

Jo

bacteria in the microbiota of different-aged workers were significantly influenced (Table 2). The relative abundance of the Gilliamella genus peaked at 36.32% at 1 dpe and then decreased to 5.91% at 3 dpe; a dynamic equilibrium was reached at 7 dpe. Snodgrassella was significantly colonized at 1 dpe and was relatively stable in the gut microbiota. The relative abundance of the genus Lactobacillus in the worker gut was only 2.79% at 0 dpe and 5.54% at 1 dpe, and the relative

abundance peaked at 48.28% at 3 dpe, indicating that the colonization of Lactobacillus occurred after three days of worker emergence. The relative abundance of Bifidobacterium in the gut of newly emerged workers was only 0.58%. The relative abundance of Bifidobacterium reached 12.10% at 3 dpe, and was stable from 19 dpe. The relative abundance of Rhizobiales was only 0.16% in the gut of newly emerged workers. After reaching 7.24% at 19 dpe, the relative abundance of Rhizobiales remained stable. The relative abundance of Frischella reached 9.70% at

ro of

1 dpe and then decreased and became stable. The relative abundance of Commensalibacter

increased significantly at 3 dpe and then also continued to increase. The relative abundance of Apibacter increased significantly at 19 dpe and then also continued to increase. The relative

-p

abundance of Bartonella increased significantly at 12 dpe. The significant increase in core

re

microbiota was accompanied by a significant decrease in genera such as Acinetobacter,

lP

Cyanobacteria, Peptostreptococcaceae, Bacteroidales, Bacteroides and Escherichia-Shigella

na

(Figure 5).

3.4 Succession of the core microbiota of workers

ur

By comparing the relative abundance of gut microbiota between the 0 dpe and 1 dpe workers, we found that the relative abundance of Acinetobacter was decreased significantly from 18.74% to

Jo

0.18% (p = 0.00508), Cyanobacteria was decreased from 6.68% to 1.50% (p = 0.00824), and Xanthomonadaceae was decreased from 1.58% to 0.015% (p = 0.00499). The relative abundance of microbiota carried by the 0 dpe workers was significantly reduced compared to the workers at other ages, and they were replaced by core microbiota. For example, Gilliamella increased from 0.0044% to 36.32% (p = 0.00367), Snodgrassella increased from 0.0045% to 9.86%

(p = 0.00716), Sphaerochaeta increased from 0.24% to 1.73% (p = 0.0202), and Campylobacter increased from 0.11% to 1.08% (p = 0.0306). In addition, Frischella increased from 0% to 9.70% (p = 0.00278). Frischella was not detected in the gut of the newly emerged workers. By comparing 1 dpe with 3 dpe, we found that this is the critical time for the colonization of Lactobacillus and Bifidobacterium because Lactobacillus increased from 5.54% to 48.28% (p = 0.00508) and Bifidobacterium increased from 0.41% to 12.10% (p = 0.00508). In addition, 1 dpe

ro of

to 3 dpe is also the key time frame of Commensalibacter colonization, which was increased from 0.021% to 2.12% (p = 0.00508). The relative abundance of the genera Gilliamella (p = 0.00824), Parabacteroides (p = 0.0202) and Alloprevotella (p = 0.0306) also decreased significantly.

-p

By comparing 7 dpe with 12 dpe, we found that Lactobacillus helsingborgensis decreased

re

from 6.08% to 1.94% (p = 0.00508). In addition, Lactobacillus kunkeei increased from 0.29% to

lP

2.20% (p = 0.00824), and Bartonella increased from 0.012% to 1.31% (p = 0.0078). By comparing 19 dpe with 25 dpe, Bifidobacterium decreased from 8.29% to 2.70% (p =

na

0.0453), and Commensalibacter increased from 10.07% to 3.32% (p = 0.0306). In addition, the three genera Bacteroides, Escherichia-Shigella, and Porphyromonadaceae, also decreased

ur

significantly during this period (Figure 6).

Jo

3.5 Beta diversity of gut bacteria Based on the weighted UniFrac distance or Bray-Curtis distance, we performed a principal coordinate analysis (PCoA) and nonmetric multidimensional scale analysis (NMDS). Additionally, an analysis of similarities (ANOSIM) was performed based on the weighted UniFrac distance. The results showed that the composition of the gut microbiota was correlated with host age. The PCoA

and NMDS plots showed that the newly emerged workers were far from the other groups, and by comparing the PCoA and NMDS plots of the other groups, the separation was also obvious at 1 dpe, indicating that the relative abundance of the core microbiota of workers was not high at 1 dpe. The 3 dpe and 7 dpe groups were closely related, and the 12 dpe and 19 dpe groups were closely related, despite the close distance of the groups at 25 dpe (group, Am25), 30 dpe, 35 dpe, and 40 dpe. Overall, samples belonging to each group were clustered in different ages of workers. In the

ro of

coordinate analysis (PCoA), PC1 accounted for 33.35% of the total variance, and PC2 accounted for 20.68%. In the nonmetric multidimensional scale analysis (NMDS), stress = 0.129, which indicated that the grouping and sampling was reliable (Figure 7). In addition, the ANOSIM

-p

analysis results were as follows: R = 0.4794, P = 0.001. This result further reveals that the

lP

re

grouping of the experiments was reasonable.

4 Discussion

na

In this study, high-throughput sequencing technology was used for deep sequencing of samples taken at multiple time (ten sample points). Our results indicate that the composition and structure

ur

of the gut flora of workers are significantly altered with host development. The core microbiome at 0 dpe was not shared with that at the other stages, which is consistent with previous studies

Jo

(Yun et al., 2018). Elijah Powell et al. found that the core microorganisms S. alvi, G. apicola, and F. perrara are significantly colonized through social contact (Powell et al., 2014). Therefore, the observation of no core gut flora in the 0 dpe worker may be explained by the minimal social contact of the workers. In addition, the relative abundance of Acinetobacter was the highest at 0 dpe, indicating that the genus Acinetobacter plays an important role in the gut of workers. The

structure of the gut flora at 0 dpe was rapidly changed at 1 dpe, particularly core genera, such as Frischella, Snodgrassella and Gilliamella, indicating that the three core genera primarily occupy the ecological niche of Acinetobacter. F. perrara causes a strong activation of the host immune system, and studies have shown that the presence of F. perrara affects gut immunity and homeostasis, and previous studies demonstrated that an increased abundance of F. perrara is related to dietary alteration and abnormal host development (Emery et al., 2017; Maes PW, 2016).

ro of

The large changes in the diet and developmental environment of the newly emerged workers probably played a main role in explaining why Frischella first colonize the worker gut. The

relative abundance of the Gilliamella genus peaked at 1 dpe, revealing a key role for Gilliamella

-p

as a core genus in the social exercise of workers. Our previous study found that the abundance of

re

G. apicola was relatively low at 1 dpe (Guo et al., 2015), which may be explained by different

lP

colonization patterns between different species of gut bacteria.

Lactobacillus and Bifidobacterium are helpful in the nutrition absorption and protection of

na

bees. The Lactobacillus and Bifidobacterium genera isolated from bee gut were sprayed into the hive, on the brood population, and on the pollen, which induced a significant promotion in

ur

harvestable honey (Alberoni et al., 2018). These two genera obviously colonized between 1 dpe and 3 dpe in this study because the relative abundance of these two species of bacteria increased

Jo

significantly during this period (from 1 dpe to 3 dpe). Anderson et al. explored Lactobacillus colonization using high-throughput sequencing technology. The colonization rule of Lactobacillus was elucidated under different diets and social contact, and the results showed that the colonization process of Lactobacillus in the hindgut fluctuated with changes in nutrition, hive and social environment (Anderson et al., 2016). Furthermore, the colonization pattern of Lactobacillus

in this study was consistent with that of previous studies (Anderson et al., 2016; Guo et al., 2015). Namely, previous studies showed that different sources did not affect the colonization mode of Lactobacillus. Between 1 dpe and 3 dpe, the abundance of the Lactobacillus and Bifidobacterium genera increased constantly. The main tasks of workers between 1 and 3 dpe include cell cleaning, resting or grooming (Seeley, 1982). Studies by Kenerova et al. (2017) showed that Bifidobacterium can stimulate the host to produce hormones, which may affect the development of

ro of

bees and accelerate the development of workers (Kešnerová et al., 2017). A significantly

increased relative abundance of Bifidobacterium may be an indicator of the rapid development of workers during this period.

-p

With the development of bee age, the social tasks of bees are diverge. Between 4 dpe and 12

re

dpe, worker bees serve as nurses and porters to transfer proteinaceous secretions to younger and

lP

older bees (Crailsheim, 1992; Crailsheim, 1991; Seeley, 1982). The relative abundance of Lactobacillus and Bifidobacterium in the gut is stable after 4 dpe. The Bifidobacterium in the bee

na

gut metabolize a range of carbohydrates, providing many nutrients and energy (Bottacini et al., 2012b) and material that guarantees social tasks for workers to perform.

ur

In addition, Lactobacillus and Bifidobacterium have many functions related to the decomposition of carbohydrates, and it has been proposed that the Lactobacillus group consisting

Jo

of Lactobacillus and Bifidobacterium is involved in the nectar processing and metabolism of carbohydrates (Butler et al., 2013; Engel et al., 2012). The colonization of these two genera promotes the digestion and absorption of nutrients in the host to indirectly promote social tasks of workers. Lactobacillus kunkeei significantly increased from 7-12 dpe, and L. kunkeei did not present a

unique distribution in the gut flora of workers. L. kunkeei has also been detected in other organs of workers. For example, plenty of L. kunkeei, an acid-resistant and osmotolerant bacterium, is found in bread (Anderson et al., 2014). Furthermore, workers at this age mainly exercised the task of honey processing (Trumbo et al., 1997). We thus speculated that L. kunkeei is a group of key bacteria in the gut microbiota of works for honey generation, and the resistance function of workers for acid and osmotolerant damage caused by honey production was formed from 7-12

ro of

dpe.

The abundance of Bifidobacterium was significantly reduced between 9 dpe and 30 dpe.

During this time, workers leave from hives and then focus on foraging for pollen, nectar, propolis

-p

and water (Calderone, 1998; Robinson, 1992). This may be one reason for the decrease in the

re

relative abundance of the Bifidobacterium genus, as reported by studies of Anderson et al. (2018),

lP

who found that the relative abundance of Bifidobacterium decreased significantly with the age of workers (Anderson et al., 2018), which is similar to results obtained in the current study. In

na

addition, F. perrara also decreases significantly due to changes in the environment and diet when workers leave the hive to conduct collection tasks (Martinson VG, 2012). This observation shows

ur

that the change in the relative abundance of F. perrara is consistent with that of Bifidobacterium during this period. However, the Bifidobacterium genus is a bacterial indicator for assessing the

Jo

age of workers (Anderson et al., 2018). The results obtained here indicate the potential of F. perrara as an indicator of age assessment, but this speculation requires further investigation in the future. In summary, as an important economic insect, the health of honey bees is of concern. High-throughput sequencing of the 16S rRNA gene was performed for the gut flora of bees across

various developmental stages in this study. We found that the structure of the gut flora of bees significantly varied with age development. The 0 dpe workers did not harbor core gut flora in the gut, and the key points for colonization of the core gut flora were around 1-3 dpe. In addition, we also elucidated the rules of dynamic microbial succession in the bee gut. The change in gut flora was significantly correlated with the increase in diurnal age. This study not only deepens our understanding of the colonization pattern of core gut flora in workers but also provides useful

ro of

information for exploring the detailed processes of bee gut flora colonization in depth.

Jo

ur

na

lP

re

-p

Acknowledgments: This study was supported by the National Natural Science Foundation of China (31660695)

Jo

ur

na

lP

re

-p

ro of

References: Alberoni, D., Baffoni, L., Gaggìa, F., Ryan, P., Murphy, K., Ross, P. R., Stanton, C., Di Gioia, D., 2018 Impact of beneficial bacteria supplementation on the gut microbiota, colony development and productivity of Apis mellifera L. Benef. Microbes 9, 269-278. Anderson, K. E., Carroll, M. J., Sheehan, T., Lanan, M. C., Mott, B. M., Maes, P., Corby-Harris, V., 2014 Hive-stored pollen of honey bees: many lines of evidence are consistent with pollen preservation, not nutrient conversion. Mol. Ecol. 23, 5904-5917. Anderson, K. E., Ricigliano, V. A., Mott, B. M., Duan, C. C., Floyd, A. S., Maes, P., 2018 The queen’s gut refines with age: longevity phenotypes in a social insect model. Microbiome 6, 108. Anderson, K. E., Rodrigues, P. A. P., Mott, B. M., Maes, P., Corby-Harris, V., 2016 Ecological Succession in the Honey Bee Gut: Shift in Lactobacillus Strain Dominance During Early Adult Development. Mol. Ecol. 71, 1008-1019. Babendreier D, J. D., Romeis J, Bigler F, Widmer F, 2007 Bacterial community structures in honeybee intestines and their response to two insecticidal proteins. FEMS Microbiol. Ecol. 59, 600 - 610. Bottacini, F., Milani, C., Turroni, F., Sã, n. B., Foroni, E., Duranti, S., Serafini, F., Viappiani, A., Strati, F., Ferrarini, A., 2012a Bifidobacterium asteroides PRL2011 genome analysis reveals clues for colonization of the insect gut. Plos One 7, e44229. Bottacini, F., Milani, C., Turroni, F., Sánchez, B., Foroni, E., Duranti, S., Serafini, F., Viappiani, A., Strati, F., Ferrarini, A., Delledonne, M., Henrissat, B., Coutinho, P., Fitzgerald, G. F., Margolles, A., van Sinderen, D., Ventura, M., 2012b Bifidobacterium asteroides PRL2011 genome analysis reveals clues for colonization of the insect gut. Plos One 7, e44229-e44229. Butler, È., Alsterfjord, M., Olofsson, T. C., Karlsson, C., Malmström, J., Vásquez, A., 2013 Proteins of novel lactic acid bacteria from Apis mellifera mellifera: an insight into the production of known extra-cellular proteins during microbial stress. BMC Microbiol. 13, 235. Calderone, N. W., 1998 Proximate mechanisms of age polyethism in the honey bee, Apis mellifera L. Apidologie 29, 127-158. Cameron, S. A., Lozier, J. D., Strange, J. P., Koch, J. B., Cordes, N., Solter, L. F., Griswold, T. L., 2011 Patterns of widespread decline in North American bumble bees. Proc. Natl. Acad. Sc.i U. S. A. 108, 662. Colin, T., Meikle, W. G., Paten, A. M., Barron, A. B., 2019 Long-term dynamics of honey bee colonies following exposure to chemical stress. Sci. Total. Environ. 677, 660-670. Crailsheim, K., 1992 The flow of jelly within a honeybee colony. J. Comp. Physiol. B 162, 681-689. Crailsheim, K., 1991 Interadult feeding of jelly in honeybee (Apis mellifera L.) colonies. J. Comp. Physiol. B 161, 55-60. Emery, O., Schmidt, K., Engel, P., 2017 Immune system stimulation by the gut symbiont Frischella perrara in the honey bee (Apis mellifera). Molecular Ecology 26, 2576-2590. Engel, P., Martinson, V. G., Moran, N. A., 2012 Functional diversity within the simple gut microbiota of the honey bee. Proc. Natl. Acad. Sc.i U. S. A. 109, 11002-11007. Guo, J., Wu, J., Chen, Y., Evans, J. D., Dai, R., Luo, W., Li, J., 2015 Characterization of gut bacteria at different developmental stages of Asian honey bees, Apis cerana. J. Invertebr. Pathol. 127, 110-114. Hauke, K., Paul, S. H., 2011 Socially transmitted gut microbiota protect bumble bees against an intestinal parasite. Proc. Natl. Acad. Sc.i U. S. A. 108, 19288-19292.

Huang SK, Y. K., Huang WF, Ying BH, Su X, Lin LH, Li JH, Chen YP, Li JL, Bao XL, Hu JZ, 2018 Influence of Feeding Type and Nosema ceranae Infection on the Gut Microbiota of Apis cerana Workers. mSystems 3, e00177-00118. Jeyaprakash, A., Hoy, M. A., Allsopp, M. H., 2003 Bacterial diversity in worker adults of Apis mellifera capensis and Apis mellifera scutellata (Insecta: Hymenoptera) assessed using 16S rRNA sequences. J. Invertebr. Pathol. 84, 96-103. Jones, J. C., Fruciano, C., Hildebrand, F., Al Toufalilia, H., Balfour, N. J., Bork, P., Engel, P., Ratnieks, F. L., Hughes, W. O., 2017 Gut microbiota composition is associated with environmental landscape in honey bees. Ecol. Evol. 8, 441-451. Kakumanu, M. L., Reeves, A. M., Anderson, T. D., Rodrigues, R. R., Williams, M. A., 2016 Honey Bee Gut Microbiome Is Altered by In-Hive Pesticide Exposures. Front. Microbiol. 7, 1255-1255. Kešnerová, L., Mars, R. A. T., Ellegaard, K. M., Troilo, M., Sauer, U., Engel, P., 2017 Disentangling metabolic functions of bacteria in the honey bee gut. PLoS Biol. 15, e2003467-e2003467.

ro of

Kwong, W. K., Moran, N. A., 2013 Cultivation and characterization of the gut symbionts of honey

bees and bumble bees: description of Snodgrassella alvi gen. nov., sp. nov., a member of the family Neisseriaceae of the Betaproteobacteria, and Gilliamella apicola gen. nov., sp. nov., a membe. International Journal of Systematic & Evolutionary Microbiology 63, 2008-2018.

Kwong, W. K., Moran, N. A., 2016 Gut microbial communities of social bees. Nat. Rev. Microbiol. 14,

-p

374.

Langowska, A., Zawilak, M., Sparks, T. H., Glazaczow, A., Tomkins, P. W., Tryjanowski, P., 2017 Long-term effect of temperature on honey yield and honeybee phenology. Int. J. Biometeorol. 61,

re

1125-1132.

Liu PD, W. J., Li HY, Lin SW, 2011 Economic values of bee pollination to China’s agriculture. Scientia Agricultura Sinica, (in chinese) 44, 5117-5123.

lP

Maes PW, R. P., Oliver R, Mott BM, Anderson KE, 2016 Diet-related gut bacterial dysbiosis correlates with impaired development, increased mortality and Nosema disease in the honeybee (Apis mellifera). Mol. Ecol. 25, 5439-5450.

Martinson, V. G., Danforth, B. N., Minckley, R. L., Olav, R., Salim, T., Moran, N. A., 2011 A simple

na

and distinctive microbiota associated with honey bees and bumble bees. Mol. Ecol. 20, 619-628. Martinson VG, M. J., Moran NA., 2012 Establishment of characteristic gut bacteria during development of the honeybee worker. Appl. Environ. Microbiol. 78, 2830-2840.

ur

Moran NA, P. G., Sandström JP, Wilcox JL, 2003 A genomic perspective on nutrient provisioning by bacterial symbionts of insects. Appl. Environ. Microbiol. 100, 14543–14548. Motta, E. V. S., Raymann, K., Moran, N. A., 2018 Glyphosate perturbs the gut microbiota of honey

Jo

bees. Proc. Natl. Acad. Sci. U. S. A. 115, 10305. Owen, R., 2017 Role of Human Action in the Spread of Honey Bee (Hymenoptera: Apidae) Pathogens. J. Econ. Entomol. 110, 797-801.

Philipp, E., Kwong, W. K., Moran, N. A., 2013 Frischella perrara gen. nov., sp. nov., a gammaproteobacterium isolated from the gut of the honeybee, Apis mellifera. Int J Syst Evol Microbiol 63, 3646-3651. Powell, J. E., Martinson, V. G., Urban-Mead, K., Moran, N. A., 2014 Routes of Acquisition of the Gut Microbiota of the Honey Bee Apis mellifera. Appl. Environ. Microb. 80, 7378-7387. Raymann K, S. Z., Moran NA., 2017 Antibiotic exposure pertubs the gut micxrobiota and elevates mortality in honeybees. PLoS Biol. 15, e2001861.

Robinson, G. E., 1992 Regulation of Division of Labor in Insect Societies. Annu. Rev. Entomol. 37, 637-665. Ross, R. P., 2012 Composition of the early intestinal microbiota:Knowledge, knowledge gaps and the use of high-throughput sequencing to address these gaps. Gut Microbes 3, 203-220. Sabree, Z. L., Hansen, A. K., Moran, N. A., 2012 Independent Studies Using Deep Sequencing Resolve the Same Set of Core Bacterial Species Dominating Gut Communities of Honey Bees. Plos One 7, e41250. Seeley, T. D., 1982 Adaptive significance of the age polyethism schedule in honeybee colonies. Behav. Ecol. Sociobio. 11, 287-293. Southwick, E. E., Southwick, L., 1992 Estimating the Economic Value of Honey Bees (Hymenoptera: Apidae) as Agricultural Pollinators in the United States. J. Econ. Entomol. 85, 621-633. Stokstad, E., 2007 The Case of the Empty Hives. Science 316, 970-972. Trumbo, S. T., Huang, Z.-Y., Robinson, G. E., 1997 Division of labor between undertaker specialists

ro of

and other middle-aged workers in honey bee colonies. Behav. Ecol. Sociobio. 41, 151-163.

Vanessa, C. H., Snyder, L. A., Schwan, M. R., Patrick, M., McFrederick, Q. S., Anderson, K. E., 2014 Origin and effect of Alpha 2.2 Acetobacteraceae in honey bee larvae and description of Parasaccharibacter apium gen. nov., sp. nov. Appl. Environ. Microbiol. 80, 7460-7472.

Yun, J. H., Jung, M. J., Kim, P. S., Bae, J. W., 2018 Social status shapes the bacterial and fungal gut

-p

communities of the honey bee. Sci. Rep. 8, 2019.

Zhang F., Sun X. X., Zhang X. C., Zhang S., Lu J., Xia Y. M., Huang Y. H., Wang X. J., 2018 The interactions between gut microbiota and entomopathogenic fungi: a potential approach for

Jo

ur

na

lP

re

biological control of Blattella germanica (L.). Pest. Manag. Sci. 74, 438-447.

Figure legends Fig:1 The figure shows the significant difference between the samples of the two groups, and

Jo

ur

na

lP

re

-p

ro of

marks the two groups with significant difference (0.01 < P ≤ 0.05 is marked as *, 0.001 < P ≤ 0.01 is marked as **, and P ≤ 0.001 is marked as ***).The abscissa is the grouping name and the ordinate is the exponential average value of each group.

Fig:2 The abscissa is the sample name, and the ordinate is the proportion of species in the sample. The column of different colors represents different species, and the length of the column represents the proportion of the species.

ro of

-p

Fig:3 The abscissa is the sample name, and the ordinate is the proportion of species in the sample.

Jo

ur

na

lP

re

The column of different colors represents different species, and the length of the column represents the proportion of the species.

Fig:4 The vertical axis represents the species names at phylum classification level, the

corresponding column length represents the average relative abundance of the species in various groups, and different colors represent different groups.On the far right is the value of P, *: 0.01 <

re

-p

ro of

P ≤ 0.05，**:0.001 < P ≤ 0.01, ***: P ≤ 0.001

lP

Fig:5 The vertical axis represents the species names at genus classification level, the corresponding column length represents the average relative abundance of the species in various groups, and different colors represent different groups.On the far right is the value of P, *: 0.01 <

Jo

ur

na

P ≤ 0.05，**:0.001 < P ≤ 0.01, ***: P ≤ 0.001

ro of -p re lP na

Fig:6Student's t-test is used to detect two different groups. The vertical axis represents the

ur

species names at genus classification level, the corresponding column length represents the average relative abundance of the species in various groups, and different colors represent

Jo

different groups.On the far right is the value of P, *: 0.01 < P ≤ 0.05，**:0.001 < P ≤ 0.01, ***: P ≤ 0.001

ro of

-p

Fig:7 The horizontal and vertical coordinates represent the two selected principal coordinate

Jo

ur

na

lP

re

components, and the percentage represents the contribution value of the principal coordinate component to the sample composition difference.The scales of the horizontal and vertical axes are relative distances and have no practical significance.Points of different colors or shapes represent samples of different groups. The closer the two sample points are, the more similar the species composition of the two samples will be.

Table 1. Richness and diversity indexes relative to each gut sample ID

Coverage

Number of

Alpha diversity

OUTs

ACE

Chao

Shannon

Simpson

0.99861

788

797

809

5.14

0.0259

0 dpe_2

0.99904

921

953

971

5.42

0.0263

0 dpe_3

0.99933

480

511

516

3.85

0.1199

0 dpe_4

0.99943

917

943

964

5.56

0.0184

0 dpe_5

0.99930

1009

1037

1052

5.53

0.0193

0 dpe_6

0.99947

782

806

817

4.69

0.0750

1 dpe_1

0.99910

398

451

458

2.46

0.2351

1 dpe_2

0.99957

510

528

532

4.12

0.0841

1 dpe_3

0.99963

557

573

576

1 dpe_4

0.99933

385

416

418

1 dpe_5

0.99949

635

670

674

1 dpe_6

0.99799

309

460

439

3 dpe_1

0.99823

154

342

315

3 dpe_2

0.99806

169

365

326

3 dpe_3

0.99704

387

647

3 dpe_4

0.99900

219

302

3 dpe_5

0.99837

264

368

3 dpe_6

0.99667

650

932

7 dpe_1

0.99781

347

526

7 dpe_2

0.99772

230

7 dpe_3

0.99791

368

7 dpe_4

0.99880

333

7 dpe_5

0.99838

150

7 dpe_6

0.99859

12 dpe_1 12 dpe_2

ro of

0 dpe_1

0.0682

2.38

0.2735

4.03

0.0778

0.83

0.7304

2.29

0.1453

2.30

0.1415

636

2.52

0.1546

307

2.59

0.1133

360

2.53

0.1194

946

3.91

0.0520

510

2.27

0.2215

436

409

2.22

0.1608

513

503

2.89

0.0942

433

439

2.51

0.1978

297

266

2.07

0.2000

264

364

374

2.58

0.1186

0.99817

184

570

322

2.11

0.1897

0.99800

250

427

405

2.55

0.1281

re

lP

na

ur

-p

4.19

0.99866

193

299

297

2.17

0.1817

12 dpe_4

0.99792

169

374

377

2.17

0.1727

12 dpe_5

0.99725

194

645

455

2.32

0.1604

Jo

12 dpe_3

12 dpe_6

0.99689

631

878

883

3.63

0.0812

19 dpe_1

0.99791

340

551

537

2.63

0.1379

19 dpe_2

0.99863

418

538

535

2.49

0.2164

19 dpe_3

0.99814

235

400

387

2.65

0.0999

19 dpe_4

0.99784

217

527

447

2.19

0.1555

19 dpe_5

0.99834

143

310

280

2.23

0.1536

19 dpe_6

0.99844

287

404

406

2.64

0.1087

25 dpe_1

0.99811

144

420

308

1.89

0.2209

25 dpe_2

0.99753

252

477

433

2.53

0.1188

0.99855

126

366

259

2.22

0.1462

25 dpe_4

0.99916

80

212

153

1.93

0.1985

25 dpe_5

0.99893

107

211

199

1.95

0.1909

25 dpe_6

0.99853

98

374

263

1.90

0.2329

30 dpe_1

0.99748

268

517

480

2.27

0.1608

30 dpe_2

0.99762

209

575

445

2.24

0.1559

30 dpe_3

0.99840

119

422

304

2.27

0.1319

30 dpe_4

0.99847

103

377

261

2.05

0.1772

30 dpe_5

0.99803

245

426

407

2.41

0.1471

30 dpe_6

0.99847

112

384

300

2.24

0.1654

35 dpe_1

0.99873

82

350

223

1.53

0.3053

35 dpe_2

0.99693

350

636

629

2.59

0.1269

35 dpe_3

0.99826

176

344

320

2.42

0.1220

35 dpe_4

0.99809

134

394

305

35 dpe_5

0.99873

99

245

252

35 dpe_6

0.99761

299

521

504

40 dpe_1

0.99792

250

446

437

40 dpe_2

0.99775

360

533

547

40 dpe_3

0.99858

104

510

255

40 dpe_4

0.99833

149

324

40 dpe_5

0.99818

167

381

40 dpe_6

0.99856

122

290

P value

-

0.000085

FDR

-

0.00043

0.1341

1.91

0.2386

2.56

0.1479

2.40

0.1440

2.74

0.1182

2.51

0.1008

297

2.53

0.1059

375

2.80

0.0781

270

2.25

0.1365

0.023

0.0024

0.0016

0.010

0.023

0.0072

0.0063

0.020

re

-p

2.28

lP

Jo

ur

na

OTUs were defined at the 97% similarity level.

ro of

25 dpe_3

f

Relative abundance (%) and time (dpe)

genera

oo

Table 2 Dynamic evolution of the top15 abundance.

P value

FDR

28.380

1.94E-05

0.000163

23.200

16.140

0.00024

0.00116

18.260

17.010

19.770

0.00846

0.0218

3.625

19.420

13.100

8.340

5.92E-05

0.000407

8.285

2.693

2.767

3.406

2.215

1.24E-05

0.000156

2.551

5.663

8.554

1.533

6.787

6.184

0.00219

0.00667

3.439

2.248

10.070

3.318

11.180

6.389

7.418

2.63E-06

5.38E-05

0.006

0.008

0.183

0.706

2.273

2.758

1.520

5.250

0.00336

0.009845

0.002

0.003

0.012

1.313

0.433

1.284

0.425

1.277

2.636

7.25E-06

0.000107

18.740

0.184

0.071

0.060

0.021

0.045

0.020

0.007

0.008

0.015

1.12E-05

0.000154

Cyanobacteria

6.677

1.499

1.257

1.257

0.904

0.906

0.138

0.148

0.115

0.044

4.81E-06

7.46E-05

Peptostreptococcaceae

1.274

3.513

2.745

0.379

2.149

0.386

0.115

0.135

0.467

0.186

0.001357

0.00436

Bacteroidales

6.381

1.975

0.460

0.298

0.370

0.247

0.043

0.119

0.118

0.159

3.94E-05

0.000295

Bacteroides

2.268

3.111

0.834

0.471

0.628

0.449

0.067

0.187

0.162

0.232

0.000146

0.000835

Escherichia-Shigella

3.621

1.917

0.994

0.207

0.661

0.677

0.046

0.058

0.124

0.097

9.53E-05

0.000624

1 dpe

3 dpe

7 dpe

12 dpe

19 dpe

25 dpe

30 dpe

35 dpe

40 dpe

Lactobacillus

2.786

5.537

48.280

40.630

36.990

30.740

26.830

27.800

24.200

Gilliamella

0.004

36.320

5.907

16.090

22.210

13.070

24.890

13.510

Snodgrassella

0.005

9.864

14.970

14.130

18.620

17.320

24.920

Rhizobiales

0.157

0.092

0.547

0.354

0.173

7.241

Bifidobacterium

0.576

0.414

12.100

13.000

5.681

Frischella

0.000

9.703

2.671

4.574

Commensalibacter

0.003

0.021

2.118

Apibacter

0

0.005

Bartonella

0

Acinetobacter

e-

Pr

na l

Jo ur

pr

0 dpe