Dynamic Changes of Gut Microbial Communities of Bumble Bee Queens through Important Life Stages.

Bumble bees are important pollinators in natural and agricultural ecosystems. Their social colonies are founded by individual queens, which, as the predominant reproductive females of colonies, contribute to colony function through worker production and fitness through male and new queen production. Therefore, queen health is paramount, but even though there has been an increasing emphasis on the role of gut microbiota for animal health, there is limited information on the gut microbial dynamics of bumble bee queens. Employing 16S rRNA amplicon sequencing and quantitative PCR, we investigate how the adult life stage and physiological state influence a queen's gut bacterial community diversity and composition in unmated, mated, and ovipositing queens of Bombus lantschouensis We found significant shifts in total gut microbe abundance and microbiota composition across queen states. There are specific compositional signatures associated with different stages, with unmated and ovipositing queens showing the greatest similarity in composition and mated queens being distinct. The bacterial genera Gilliamella, Snodgrassella, and Lactobacillus were relatively dominant in unmated and ovipositing queens, with Bifidobacterium dominant in ovipositing queens only. Bacillus, Lactococcus, and Pseudomonas increased following queen mating. Intriguingly, however, further analysis of unmated queens matching the mated queens in age showed that changes are independent of the act of mating. Our study is the first to explore the gut microbiome of bumble bee queens across key life stages from adult eclosion to egg laying and provides useful information for future studies of the function of gut bacteria in queen development and colony performance.IMPORTANCE Bumble bee queens undergo a number of biological changes as they transition through adult emergence, mating, overwintering, foraging, and colony initiation including egg laying. Therefore, they represent an important system to understand the link between physiological, behavioral, and environmental changes and host-associated microbiota. It is plausible that the bumble bee queen gut bacteria play a role in shaping the ability of the queen to survive environmental extremes and reproduce, due to long-established coevolutionary relationships between the host and microbiome members.

INTRODUCTION

Our results show that there is a significant difference in diversity and composition of the gut microbial communities in queens of Bombus lantschouensis across different states. This study provides insight into the relationship between the bacterial community and the physiological states of bumble bee queens and lays the foundation for further studies of the functioning of the gut microbiota in the health and reproductive success of bumble bee queens.

Symbiotic bacteria play important roles in physiology, behavior, and ultimately fitness of their animal hosts, including insects (1, 2). They can supply the necessary nutrition for their host (3, 4), improve the host’s development and fecundity (5, 6), modulate their metabolism (7), and induce insects’ aggregations (8, 9) as well as determine kin recognition and mate choice (10–12). Conversely, host physiology and developmental stage can influence the composition of host-associated microbiota (13, 14). Bumble bees are social insects with annual colony cycles, and queens undergo a number of physiological and developmental changes as they progress through mating and diapause to subsequent oviposition and colony production (15, 16). However, there is limited information on the association between these changes and the critical gut microbiota of bumble bees.

Compared with the gut microbiota of many other animals, adult workers of social apid bees (bumble bees and honey bees) harbor a relatively simple yet specialized gut microbiota dominated by several recently described bacterial species, including Gilliamella apicola, Snodgrassella alvi, and specialized species of Lactobacillus (17–22). A number of beneficial interactions among these microbes and the honey and bumble bees have been reported, including increased metabolic functionality, protection from invading pathogens through facilitation of the immune response (23, 24), or exclusionary effects (25–27).

In contrast to worker bees, few studies have examined microbial communities that are associated with honey bee and bumble bee queens, even though their health and proper function are central to the productivity of their colonies. Parmentier et al. (28) found that typical core gut microbial communities in adults are absent in the larvae of wild bumble bees, which suggested that the compositions of microbial communities are different among different developmental stages or castes of bumble bees. The microbiota has also been shown to change during the developmental trajectory of honey bee queens (29) and the hibernation of bumble bee queens (30). Tarpy et al. (29) suggested that mainly enteric bacteria are present in honey bee queens at an early stage, with the predominant gut bacteria being Alphaproteobacteria at maturity, and yet the size and composition of workers’ symbiotic bacteria are relatively stable across ages. Likewise, through an isolated queen experiment, Powell et al. (31) found fewer of the primary honey bee microbiota in honey bee queen guts, and the species of bacteria were different from those of workers. Queens of honey bees and bumble bees have distinct biology. The latter goes through diapause and has a solitary founding stage, neither of which occurs in honey bees. These differences and the fact that founding bumble bee queens are the source of certain components of the gut microbiota of workers in the subsequent colony (32) make understanding the dynamics of the gut microbiota across key life stages of bumble bee queens crucial. In bumble bee queens, Bosmans et al. (30) revealed that the bacterial community composition during hibernation is richer, including some psychrophilic and psychrotrophic taxa, than in nonhibernating active queens. These studies indicate that changes of gut microbiota of honey bee and bumble bee queens may be associated with the physiological variation and developmental stage. However, temporal dynamics of bumble bee queen gut microbiota remain underexplored during sexual maturity, and patterns may provide novel insights into the interplay between queen development and physiology and the queen’s microbiota, the change to which may offer feedback on microbiota functioning in queen hosts.

Bumble bees are important pollinators of many flowering plant species in temperate to subarctic and alpine zones (33). In recent decades, many bumble bee species have been identified as declining, particularly in Europe and North America (34). Many factors have been suggested to be responsible for these declines (35). With ongoing land use and climate changes, some bumble bee species have also been predicted to become critically endangered and vulnerable in central mainland China and northeastern Asia in the future (36). At the same time, techniques for artificial hibernation and large-scale propagation have been developed that enable the commercial production of bumble bee colonies in the hundreds of thousands annually (37). Population declines and the agricultural importance of bumble bees necessitate a greater understanding of factors implicated in bumble bee health, such as the gut microbiota and their health-related functions.

Temperate bumble bees have an annual eusocial life cycle, with a solitary queen phase between mating and the foundation of new colonies in spring following the emergence of queens from hibernation (16). Toward the end of the colony cycle in late summer, sexuals (virgin queens and males) are produced. Young queens mate, usually with a single male for many species; hibernate; and subsequently emerge to produce the next generation (Fig. 1) (38, 39). The queen is critical to the development of the microbiota of individuals within the colony, with vertical transmission of certain core gut microbes occurring from mother to offspring (32). The composition of the queen’s microbiota and its contribution to health are also critically important given that the queen is the principal reproductive female in a colony, playing a crucial role in colony development and longevity (40–42).

FIG 1

Bumble bee life cycle. In the wild, Bombus lantschouensis queens emerge from hibernation in spring, forage (A), and find a nesting location in which to lay eggs and initiate a new colony (B). For the initiation of colonies to produce experimental queens in the laboratory, spring queens were collected from the field (first generation). The worker population grows, and toward the end of the colony cycle in late summer, sexuals (virgin queens and males) are produced (C). Young queens mate with only one male (D) and subsequently hibernate to produce the next generation (E, second generation B). In this study, we assessed three stages of Bombus lantschouensis queens: unmated queens (UQs; virgin queens) following queen eclosion to adulthood in colonies (C), mated queens (MQs; mating successfully with drone, D), and ovipositing queens (OQs; queens actively laying eggs, second generation B).

A better understanding of the composition of gut microbiota in different physiological states of bumble bee queens would shed light on the complex interplay between the microbiota and queen health. To this end, using an amplicon sequencing approach, we assessed microbiota composition in three queen types of the bumble bee Bombus lantschouensis: unmated queens (UQs; virgin queens) 7 days after adult eclosion; mated queens (MQs), mated at 7 days posteclosion and sampled 7 days later; and posthibernation ovipositing queens (OQs; queens actively laying eggs). These queens were offspring from laboratory-reared colonies (Fig. 1) and fed a controlled diet. To confirm differences, targeted quantification PCR (qPCR) was carried out on dominant bacterial genera identified in the three groups. This analysis included a temporal analysis of abundances comparing unmated (1 to 7 days posteclosion) and mated (1 to 7 days postmating) queens. To further distinguish intrinsic changes from those associated with mating itself, abundances were also assessed in another group of unmated queens sampled every 2 days between 1 and 15 days posteclosion. Our study is the first to explore dynamic changes of the gut microbiota across important life stages of bumble bee queens, from adult eclosion to colony foundation and egg laying. It shows dynamic diversity and variation of gut bacterial communities and improves our understanding of possible relationships between the gut microbial communities and different developmental and physiological states of bumble bee queens.

RESULTS

16S rRNA gene sequencing analysis and taxa generated.

We processed and filtered sequences, clustered them into operational taxonomic units (OTUs) with 97% minimum identity, and excluded plastids, singletons, and OTUs restricted to single samples. A total of 2,107,642 sequences were obtained in 86 samples including unmated queens (UQs; n = 30), mated queens (MQs; n = 27), and ovipositing queens (OQs; n = 29), which range from 9,915 to 44,451 (24,507 ± 961, mean ± standard error [SE]) per sample. These sequences were clustered into 390 OTUs, with a range of 17 to 138 per sample (67 ± 4, mean ± SE). The core OTUs comprised approximately 17.94% of the total candidates, while 282, 186, and 175 OTUs were identified uniquely in the UQ, MQ, and OQ groups, respectively (see Fig. S1 in the supplemental material).

Common and specific OTU distribution of gut microbiota among three groups of queens. Download FIG S1, PDF file, 0.05 MB.

To better understand the differences in the microbiome between the three queen stages, OTU sequences were blasted against the annotated SILVA 16S rRNA reference database (https://www.arb-silva.de). Twenty phyla were detected across all samples. However, four bacterial phyla accounted for more than 99% of all sequences. Ranked by relative abundance, these phyla were Proteobacteria (66%), Firmicutes (26%), Actinobacteria (6%), and Bacteroidetes (1%). Unclassified bacteria at the phylum level were rare and represented less than 1% of all sequences (Fig. S2).

Relative abundance of predominant gut bacterial phyla in three groups of queens. Download FIG S2, PDF file, 0.03 MB.

Estimates and variations of microbial local diversity among samples from three stages of queens.

We employed three species-richness measures of richness, Chao1, and abundance-based coverage estimator (ACE) to investigate the number of different OTUs (i.e., species richness) between queen groups. All measures gave qualitatively similar results, with the richness significantly affected by group identity (Fig. 2). Gut microbial communities from the MQ state had the highest richness, followed by UQs, and OQs with the lowest (P < 0.0001). Evenness was also calculated for the microbial communities, to investigate how the queen status influenced the equality of distribution of different microbes within each gut. The highest species evenness was observed for the MQ state, indicating that abundances of the diverse gut microbes associated with the MQ state were the most evenly distributed (Fig. 2). Unlike richness, the evenness of communities of OQs was greater than those of UQs. Finally, Simpson and Shannon diversity indices were calculated incorporating richness and evenness. The results for overall diversity mirrored those for evenness, with MQs having the greatest diversity, followed by OQs, and UQs having the lowest (Fig. 2). The evenness and diversity results indicate that the gut microbial community of the unmated queen is dominated by only a few species. Indeed, the bacterial genera Gilliamella and Snodgrassella were the two most dominant gut microbiota taxa, accounting for 82.5% ± 3.49% (mean ± SE) of all total sequence reads in unmated queens. In brief, our results suggest variation in the gut microbial community structure of the three queen stages (UQ, MQ, and OQ).

FIG 2

The microbial community diversity of different queen states. Box plots show OTU measures of raw richness, ACE, Chao1, evenness, and Simpson and Shannon diversity indices. *** indicates significant differences among groups (P < 0.0001). UQ, unmated queen; MQ, mated queen; OQ, ovipositing queen.

To test the above inference, we conducted beta-diversity analysis of the microbial communities among UQ, MQ, and OQ states using both unsupervised and supervised methods. The unsupervised nonmetric multidimensional scaling (NMDS) plots (Bray-Curtis distance matrix) revealed a clear separation of samples according to the different bumble bee queen stages (stress = 0.101 provides a good representation in reduced dimensions) (Fig. 3A). The supervised redundancy analysis (RDA) further indicated that as the queens progress through the different physiological states, these can significantly influence the gut microbial community composition of the host (P = 0.001, Fig. 3B). These results obtained by two independent methods to assess beta-diversity consistently suggested that there were significant structural separations of the gut microbiota among UQs, MQs, and OQs. The microbial community structure associated with the MQ state was significantly distinguished from UQ and OQ states (implicated by RDA1 with 55% of variation in Fig. 3B). Enterotype analysis of the genus-level table for the microbial communities of all samples also resulted in an optimal number of enterotypes (clusters of similar communities) of three (Fig. S3). This finding indicates that the gut microbial community structure of mated queens is unique and may be a consequence of physiological changes associated with mating itself or the development of the microbiota as the queen ages or moves toward diapause. However, the greater similarity between UQs and OQs suggests that the shift is more likely to be associated with changes during the time or process of mating.

FIG 3

The similarity and variation of gut microbial community structures across the three queen groups of unmated queens (UQ), mated queens (MQ), and ovipositing queens (OQ). (A) Unsupervised NMDS plot of beta-diversity (Bray-Curtis dissimilarity) showing clustering of the gut microbiota from individual samples (n = 86). Distances between individual samples reflect the extent of the similarity of gut microbiota. (B) Supervised RDA of the relationship between queen states and the relative abundance of bacterial genera. The P value above the figure indicates that the variable (UQ, MQ, and OQ) significantly explains variation in sample distances.

Principal-component analysis (PCA) plot of enterotype clusters observed in Bombus lantschouensis on three groups of queens (n = 86). Download FIG S3, PDF file, 0.4 MB.

The discovery of gut microbial biomarkers associated with different queen states.

Specific gut microbiota features can reflect specific disease or normal physiological conditions of the host (43). Hence, we analyzed the relative abundance of four phyla and their distribution in the guts of three stages of queens (Fig. S2). Proteobacteria ware the most predominant gut bacterial phylum in both UQs (96.4% ± 0.77%, mean ± SE) and OQs (72.2% ± 4.42%, mean ± SE). In contrast, Firmicutes account for 65.5% ± 1.16% (mean ± SE) of total reads in MQs. These results indicate potential broad-scale gut microbial markers unique to the three queen stages. Using LEfSe analysis, we identified gut microbiota features specific to the three queen groups (Fig. 4A). Among them, there were four bacterial genera in UQs, five genera in OQs, and 23 genera in MQs (Fig. 4A). Combined with a heatmap view of the relative abundance of gut bacterial genera, our results show that seven bacterial genera can be used as signatures of the three queen states (Fig. 4B). These genera showed significant differences in relative abundance between queen states (P < 0.05, Kruskal-Wallis test). Gilliamella and Snodgrassella were abundant bacterial genera associated with both UQs and OQs, with high abundances of Lactobacillus and Bifidobacterium also associated with OQs. The main bacterial genera in mated queen (MQ) guts were Bacillus, Lactococcus, and Pseudomonas. Besides the identified highly abundant bacterial phylotypes for the MQ state, there were also low-abundance bacterial genera associated with MQs, including Proteobacteria (Psychrobacter, Methylobacterium, Serratia, Escherichia, Comamonas, Citrobacter, Janthinobacterium, Stenotrophomonas, Brevundimonas, Hafnia, and Acinetobacter), Firmicutes (Brochothrix, Oceanobacillus, Geobacillus, Solibacillus, Lysinibacillus, Carnobacterium, Enterococcus, Streptococcus, and Clostridium sensu stricto), Actinobacteria (Arthrobacter and Paeniglutamicibacter), and Bacteroidetes (Myroides, Chryseobacterium, and Flavobacterium) (Fig. 4). The relatively diverse compositions in MQs were consistent with the finding of the highest alpha-diversity being in MQs, as presented in Fig. 2. These results demonstrate specific gut microbial features associated with queens at the postmating stage.

FIG 4

Comparison of the predominant bacterial genera among the three queen groupings of unmated queens (UQ), mated queens (MQ), and ovipositing queens (OQ). (A) LEfSe analysis indicates significantly different abundances of bacterial genera for each group. LDA score value is 2. (B) Heatmap depicts the relative abundances of the identified gut bacterial genera across the different queen groupings. Relative abundance of the gut bacterial genera Gilliamella, Snodgrassella, Lactobacillus, Bifidobacterium, Bacillus, Pseudomonas, and Lactococcus among three groups with significant differences (P < 0.05, Kruskal-Wallis tests). qPCR analyses were used to validate relative abundances of the seven differentially abundant genera marked by red arrowheads.

Copy number validation of differentially abundant bacterial genera in the three queen stages.

To verify the accuracy of culture-independent analysis of bacterial genera described by Fig. 4, we used 16S rRNA gene-targeted group-specific primers for real-time PCR analysis of seven identified predominant bacterial genera (Gilliamella, Snodgrassella, Lactobacillus, Bifidobacterium, Bacillus, Lactococcus, and Pseudomonas) in unmated queens, mated queens, and ovipositing queens. The mean absolute number (± SD) of the overall bacterial rRNA genes of each queen stage and age was estimated to vary from 1.47 × 108 (± [4.5 × 107]) to 6.35 × 108 (± [1.1 × 108]) copies per gut, with each of the predominant bacterial genera differentially contributing to this total (Fig. 5 and 6).

FIG 5

The change of copy numbers of overall bacteria (16S rRNA copies) and bacterial genera at different time points of the three queen groupings. UQ, unmated queen; MQ, mated queens; OQ, ovipositing queens. Days (d) for UQ represent days posteclosion and for MQ represent days postmating, with mating occurring at 7 days posteclosion. Bars represent means ± SEM. Different letters above bars within plots represent significant differences (pairwise t tests, P < 0.05).

FIG 6

The change of copy numbers of overall bacteria (16S rRNA copies) and bacterial genera at different time points for UQs (unmated queens). Days (d) represent days posteclosion. Bars represent means ± SEM. Different letters above bars within plots represent significant differences (pairwise t tests, P < 0.05).

This targeted approach confirmed the patterns seen in Fig. 4 for UQs 7 days posteclosion, MQs 7 days after mating, and OQs (Fig. 5). In addition, we investigated the temporal dynamics of changes within UQ, OQ, and MQ groups. Gilliamella, Snodgrassella, and Lactobacillus increased with age posteclosion (1 to 7 days) in UQs and declined significantly following mating at 7 days across the MQ stage (1 to 7 days postmating) but rebounded to the peak seen in UQs 7 days posteclosion by the time of queen egg laying (OQ). This initial increase, decline, and rebound were also present in the total 16S rRNA copies (Fig. 5), indicating changes in total bacterial abundance. However, the decline from 7 days posteclosion was not as pronounced due to increases in other genera. Bacillus, Lactococcus, and Pseudomonas increased postmating from low levels in the UQ state, returning to UQ state abundances in OQs. Uniquely, the bacterial genus Bifidobacterium was found at higher abundances only in OQs.

To further understand the relationship of microbiota composition with queen state and age, the same bacterial genera and total bacterial abundance were quantified in unmated queens at eight time points from 1 to 15 days posteclosion (Fig. 6). Interestingly, the results for all genera qualitatively mirrored those across UQ and MQ states shown in Fig. 5, even though these queens remained unmated throughout the time period. This indicates that the changes in microbiota composition in mated queens are not entirely the result of mating per se but rather are the consequence of age- or physiology-related changes as queens reach maturity and approach hibernation.

DISCUSSION

Our study demonstrates temporal dynamics in the gut microbiota of bumble bee queens at 1 to 15 days post-adult eclosion and posthibernation, comparing three physiological states: unmated, mated, and ovipositing. Over the first week following eclosion, we see increases in three well-known apid bee symbionts, Gilliamella, Snodgrassella, and Lactobacillus, which have previously been reported in worker bumble bees (44) and also show such a similar temporal increase in honey bee workers (45). The prehibernation peak of these three bacterial genera is consistent with sexual maturity in bumble bee queens, indicating their potential to functionally contribute to the promotion of physiological development of virgin bumble bee queens. Surprisingly, independent of mating status, after the 7-day peak, these early-abundant three bacterial genera were dramatically replaced by other bacteria, such as those members belonging to the Bacillus genus. The independence of this change from actual mating suggests that the shift may be the consequence of age- or physiology-related changes as queens reach maturity and approach hibernation. This change is particularly interesting given that the predominance of Gilliamella, Snodgrassella, and Lactobacillus was restored in posthibernation ovipositing queens. Additionally, an increase in Bifidobacterium is unique to the ovipositing queen group. Bifidobacterium increases during late pregnancy in humans, with a potential beneficial role (46, 47). Overall, the dynamic nature of the microbiota of bumble bee queens, including of some core and functionally important taxa, suggests links to the specific biology of bumble bee queens.

Whether gut microbiotas influence sexual maturation of their animal hosts remains little explored. Our results imply a connection between a change in the microbiota of bumble bee queens and the period of sexual maturity, but causation cannot be inferred. It has been shown that the mouse microbiota is required for sex-specific diurnal rhythms of gene expression and metabolism, showing that it plays a key role in ensuring proper sexual maturation and growth hormone secretion (48). In addition, studies have found that gut symbionts have potential effects on reproductive behaviors in insects. For example, in Drosophila melanogaster commensal bacteria play a role in mating preferences (10) and alteration of female microbiota counteracts a default male outbreeding strategy by inhibiting female sexual signaling (49). The gut microbiota has also been shown to modify olfactory sense-guided microbial preferences and foraging decisions in Drosophila, indicating a role of animal microbiota in shaping host fitness-related behavior through their chemosensory responses (50). Moreover, the gut bacterium Lactobacillus brevis has been found to modulate locomotor activity in D. melanogaster, which is mediated by the level of a sugar and the activity of neurons that produce the molecule octopamine (51). While the changes in the microbiota around the time of mating that we have uncovered in our study are intriguing, further work would be required to elucidate if the microbiota influences bumble bee queen mating behavior and chemical communication required for copulation, such as queen sex pheromone production (52).

Disruption of the gut microbiota of primary termite reproductives has been shown to have negative consequences for reproduction (53). A similar connection has been made between microbiota presence and parthenogenetic reproduction in Daphnia water fleas (54). Unlike these experimentally induced disturbances of the microbiota, our results show a significant but naturally occurring change in the bumble bee queen microbiota. State-related changes have been shown for humans, with a shift in the gut microbiota during pregnancy (55–57), but our study is one of the first to report such a dramatic shift in an insect. These alterations could be adaptive, with positive effects on physiological development and behavior, or simply a side effect of intrinsic physiological changes ongoing within queens as they mature or of interactions with diet (58). Of importance for understanding the causes and consequences of these dynamic changes is the demonstration that while changes occurred around the time of mating, they are not the result of mating itself. However, further investigations are required to infer causation and any consequences of the decrease of earlier core bacteria and enrichment of Bacillus, Pseudomonas, and Lactococcus in the mature bumble bee queen. There is a potential for a direct active role of these enriched genera. Sabaté et al. (59) showed that Bacillus subtilis strains isolated from honey bee gut produced surfactins and fungicides that can inhibit the important honey bee pathogens Paenibacillus larvae and Ascosphaera apis. Also, Bacillus species can produce amylase that helps in the processing of flower nectar into honey in honey bees (60). In rotifers, Lactococcus was found to serve as a probiotic to enhance growth and immunity (61, 62). Pseudomonas species in insects have been shown to be involved in detoxification (63) and digestion through amylolytic, xylanolytic, and diazotrophic activities that could contribute to the nutritional supplement and nitrogen balance (64, 65).

A particularly intriguing finding is that the changes seen in the microbiota after 7 days posteclosion in both mated and unmated queens are reversed in ovipositing queens posthibernation. This restoration of the dominance of core gut microbes including Gilliamella, Snodgrassella, and Lactobacillus, in addition to an increase in Bifidobacterium, is likely critical for colony success, given the key role that many of these taxa play in bees. Hibernation constitutes a period of considerable environmental and physiological changes, yet relatively little is known about the relationship between gut microbiota and hibernation. However, hibernation has been shown to be associated with changes in the microbiota of some organisms (66–68). For example, Sommer et al. (68) showed that the microbiota and serum metabolites in brown bears differ seasonally between hibernating and active phases and that transplants of the specific microbiota into mice transferred some of the seasonal metabolic features seen in bears. For temperate bumble bees, hibernation is usual after mating and before oviposition as an adaptation to challenges imposed by winter. This period of diapause is associated with many changes in metabolism and physiology in general, alongside the environmental alterations (15). A difference in the microbiota between queens before and during hibernation has been observed (30), which could be associated with the period of hibernation itself. However, our observations indicate that significant changes to the gut microbial community of bumble bee queens occur prior to entrance into hibernation and are largely reset in queens following hibernation, when they are egg laying. Bumble bee queens utilize storage of energy, such as lipids and glycogen, to survive low winter temperatures (69). Similar to the work of Bosmans et al. (30), we detected some cold-loving and cold-tolerant bacterial genera, such as Acinetobacter, Chryseobacterium, Hafnia, Psychrobacter, and Pseudomonas, in samples of mated queens and also older unmated queens prior to hibernation. Their presence, even prior to the initiation of abiotic environmental changes associated with hibernation, could support the queens during the environmental transition during hibernation. The contribution of this distinct microbial community to hibernation success, relative to earlier and later microbial community compositions, is thus important to investigate further.

The microbiotas associated with organisms may be closely linked with physiological and behavioral changes in their hosts, either responding to these changes indirectly or directly being involved in them. Bumble bee queens undergo a number of biological changes as they transition through adult emergence, mating, overwintering, foraging, and colony initiation including egg laying. Therefore, they represent an important system to understand the link between physiological, behavioral, and environmental changes and host-associated microbiota. It is plausible that the bumble bee queen gut bacteria play a role in shaping the ability of the queen to survive environmental extremes and reproduce, due to long-established coevolutionary relationships between the host and microbiome members. Our results show that there is a significant difference in diversity and composition of the gut bacterial species between bumble bee queens at different ages and physiological stages. This provides a critical insight into the relationship between the bacterial community and queen status in bumble bees and establishes the basis for further work to determine if the microbiota changes identified are causal in the health and success of critically important bumble bee queens.

MATERIALS AND METHODS

Overview of sampling.

An approach using Illumina amplicon sequencing of the V3-V4 region of bacterial 16S rRNA was used to investigate differences in the microbiomes of queens at the different developmental stages of virgin unmated queens (UQs; 7 days post-adult eclosion, n = 30), mated queens (MQs; mated at 7 days posteclosion and sampled 7 days later, n = 27), and ovipositing queens after diapause (OQs; n = 29). A targeted approach focusing on bacterial genera to verify these results by qPCR and to additionally assess temporal changes was carried out in UQs at 1, 3, 5, and 7 days post-adult eclosion and in MQs at 1, 3, 5, and 7 days postmating, taking place 7 days posteclosion (n = 5 per time point). Furthermore, to investigate changes independent of mating, a final temporal assessment of UQs was carried out covering the period when MQs were sampled (1, 3, 5, 7, 9, 11, 13, and 15 days posteclosion, n = 5 per time point).

Queen stages and sample collection.

Queens of B. lantschouensis were collected in the spring of 2016 from natural populations in Gansu and Ningxia provinces of China and identified by morphology and molecular methods (70, 71). Since the animals investigated in this study are neither vertebrates nor regulated invertebrates, ethical approval was not required, and the bees were sampled on property in Gansu and Ningxia provinces with consent of the manager of the Botanical Garden. Collected queens were reared in small plastic cages in the dark at a temperature of 27 ± 1°C and relative humidity of 50 to 60%. Sugar water (1:1, vol/vol) and apricot pollen were provided ad libitum to 100 colonies subsequently produced until males and gynes (new queens) emerged. Sampling for UQs, MQs, and OQs was carried out as outlined above. For matings, at 7 days post-queen eclosion, queens and males were kept at a ratio of 1:2, respectively, in a 4-m by 3-m by 2-m (length by width by height) net enclosure to ensure that one queen would mate with one male. For sampling of postdiapause OQs, mated queens were reared in a small wooden box until they became less active; they were then transferred to 4°C for diapause. After 4 months, they were revived and fed in the dark under the rearing conditions described above. Queens that had laid eggs and whose first batch of workers had emerged were collected for the OQ group. All collected samples were snap-frozen in liquid nitrogen and then stored at −80°C until the subsequent molecular analyses.

Extraction of the gut DNA.

Before removing the whole gut from queens, including crop, midgut, ileum, and rectum, the sample surface was sterilized individually with 70% and 90% ethanol solution for 1 min, respectively, followed by multiple washes using double-distilled water. The abdomen was dissected with sterilized scissors and tweezers, and the whole alimentary canal was removed and transferred into a 1.5-ml microcentrifuge tube filled with 100 μl double-distilled water and ceramic beads (0.1 mm) for the subsequent DNA extraction.

Gut samples were homogenized in a tissue lyser (Qiagen, Hilden, Germany) followed by genomic DNA isolation using the Wizard Genomic DNA purification kit (Promega; A1120) according to the manufacturer’s instructions, with DNA suspended in 30 μl nuclease-free water. The concentration and quality of extracted DNA were assessed using a Qubit fluorometer (Invitrogen, Carlsbad, CA, USA) and 2% agarose gel electrophoresis, respectively. Extracted DNA was stored at −20°C until further processing.

Illumina sequencing and bioinformatics analysis.

The hypervariable V3-V4 region of the bacterial 16S rRNA gene was amplified with the primers 341F (5′-CCTAYGGGRBGCASCAG-3′) and 806R (5′- GGACTACNNGGGTATCTAAT-3′) (72). Twenty-microliter PCR mixtures were set up with 4 μl 5× FastPfu buffer, 2 μl deoxynucleoside triphosphates (dNTPs) (2.5 mM), 0.8 μl each primer, 0.4 μl FastPfu polymerase, and template DNA (10 ng). Reactions proceeded in a GeneAmp 9700 (ABI) thermocycler with 95°C for 5 min; 27 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and elongation at 72°C for 45 s, followed by an additional elongation at 72°C for 10 min; and a dissociation stage at the end of the run. PCR products were detected by 2% agarose gel electrophoresis and purified using the QIAquick gel extraction kit (Qiagen). Library pools were constructed with equal amounts of each PCR product by using the TruSeq Nano DNA LT sample prep kit (Illumina), which was amplified through paired-end sequence on the Illumina MiSeq platform.

Raw Illumina sequence reads were modified by filtration, merging, and quality control, and barcode and primer sequences were removed, leaving library-specific tags. The fastq-join method was used to merge sequences using QIIME software (73), with an overlap length larger than 10 bp and mismatch ratio lower than 20%. Operational taxonomic unit (OTU) analysis was performed using the Uparse package (version 7.0.1001) with a 97% sequence identity on the basis of the effective tags (74). Each OTU was taxonomically assigned based on the SILVA 16S rRNA reference database using the assign_taxonomy.py program (http://qiime.org/scripts/assign_taxonomy.html). OTUs were processed by removing chloroplast sequences, mitochondrial sequences, and unclassified sequences and then obtaining species annotation information (confident threshold value, >0.8) (75, 76). Statistical differences in relative abundances of OTUs in different samples were analyzed by a nonparametric Kruskal-Wallis test, with analyses carried out using SPSS (version 17). The OTUs with relative abundance values of >0.001% (above three tags per sample) in at least one sample were retained.

The bacterial community diversities of gut samples were calculated and analyzed using the online software Calypso (http://cgenome.net/wiki/index.php/Calypso) with square root-based normalization of relative abundance. After samples were rarefied to even read depth, alpha-diversity measures of richness (Chao1 and ACE), evenness, and diversity indices (Simpson and Shannon) were compared between unmated, mated, and ovipositing queens by ANOVA. To determine if there are significant difference of gut microbial community structures among the three bumble bee queen states, the unsupervised nonmetric multidimensional scaling (NMDS) analysis of beta-diversity (Bray-Curtis distance matrix) was first conducted (77). The supervised redundancy analysis (RDA) was used to further validate complex associations between community composition and multiple explanatory variables (i.e., unmated, mated, and ovipositing queens in our study) (78). The P value reported indicates if each explanatory variable is significantly associated with variation in gut microbial composition. Gut bacterial genera associated with different physiological conditions of bumble bee queens were further identified using the linear discriminant analysis (LDA) effect size method (LEfSe) with default parameters (79). Enterotype analysis was carried out as in a previous study (44).

Genus-specific primer design and PCR amplification.

16S rRNA gene sequences of key bacterial genera were retrieved from the GenBank database. The software DNAMAN was used to align and analyze sequences and identify highly conserved regions for designing primer pairs that were unique for each genus, using Primer Premier, version 5.0. Primer sequences of Bacillus, Pseudomonas, and Lactococcus were BacF (GATGCGTAGCCGACCTGAGA) and BacR (GGCGTTGCTCCGTCAGACTT), PseF (CCGTAACTGGTCTGAGAGGATG) and PseR (GCATGGCTGGATCAGGCTTT), and LactF (GCGATGATACATAGCCGACCTG) and LactR (AGTTAGCCGTCCCTTTCTGGTT), respectively. Primers for Gilliamella, Snodgrassella, Lactobacillus, and Bifidobacterium were from the previous study (80–82), as were the universal 16S rRNA primers to determine overall bacterial load in each queen gut sample (83, 84). To confirm the specificity of each bacterial primer set, 20-μl PCR amplification was performed in a reaction mixture containing 10 μl of SYBR Premix Ex Taq II (Tli RNase H Plus) (2×), 0.8 μl of the forward primer (10 μM), 0.8 μl of the reverse primer (10 μM), 1 μl of DNA sample, and 7.4 μl of double-distilled water. The PCR cycling conditions were as follows: predenaturation at 95°C for 30 s, followed by 40 cycles of denaturation at 95°C for 5 s and annealing at 60°C for 30 s, with a subsequent melt curve to check the specificity of the amplified fragments. The product sizes of PCR amplification were confirmed by 2% agarose gel electrophoresis.

Absolute qPCR assay.

Single-band PCR products were purified using the EasyPure PCR purification kit and then inserted into the T vector using the pEASY-T1 simple cloning kit. The recombinant plasmid DNA was transformed into competent cells. Mixtures were uniformly smeared on Luria broth (LB) agar plates and cultured overnight at 37°C. The positive bacterial clones were selected and used for plasmid extraction according to the manufacturer’s instructions for the AxyPrep plasmid DNA minikit (Axygen; APMNP50). The concentration and quality of recombinant plasmid were measured by spectrophotometry (NanoDrop 2000, ThermoFisher) and visualized through 2% agarose gel electrophoresis, respectively. The recombinant plasmid DNA was stored at −80°C until use.

Based on the concentration of recombinant plasmids and the formula developed by Dhanasekaran et al. (85), copy numbers of the recombinant plasmid DNA were determined. The stock solution was 10-fold serially diluted to achieve different concentrations (from 108 to 103 copies/μl) to generate a standard curve.

Absolute quantitative PCR was performed in parallel with samples and corresponding serially diluted standards. The reaction mixture and thermocycler conditions of the PCR were the same as described above. Template DNA was diluted 10-fold before use. Each sample was run in triplicate. The actual copy numbers of specific bacterial 16S rRNA genes in samples were calculated by the threshold cycle (C) value relative to the relevant standard curve (86). Each standard curve was constructed by a liner regression of the logarithmic values of the estimated copy number of diluted standards (x axis) against the corresponding C values (y axis). The amplification efficiency (E) was related to the slope according to the formula E = 10(−1/slope) − 1 (87). The analyses of genus-specific bacterial 16S rRNA gene copy numbers in different samples were performed in SPSS (version 17). Values were normalized with log transformation. The significant differences in the copy numbers of bacteria at different time points were determined by one-way ANOVAs and least significant difference tests on the log-transformed values.

Data availability.

The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (88) in the BIG Data Center (89), Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under accession number CRA001462 and are publicly accessible at http://bigd.big.ac.cn/gsa.