The Queen Gut Refines with Age: Longevity Phenotypes in a Social Insect Model

Background In social insects, identical genotypes can show extreme lifespan variation providing a unique perspective on age-associated microbial succession. In honey bees, short and long-lived host phenotypes are polarized by a suite of age-associated factors including hormones, nutrition, immune senescence and oxidative stress. Similar to other model organisms, the aging gut microbiota of short-lived (worker) honey bees accrue Proteobacteria and are depleted of Lactobacillus and Bifidobacterium, consistent with a suite of host senescence markers. In contrast, long-lived (queen) honey bees maintain youthful cellular function without expressing oxidative stress genes, suggesting a very different host environment for age-associated microbial succession. Results We sequenced the microbiota of 63 honey bee queens exploring two chronological ages and four alimentary tract niches. To control for individual variation we quantified carbonyl accumulation in queen fat body tissue as a proxy for biological aging. We compared our results to the age-specific microbial succession of worker guts. Accounting for queen source variation, two or more bacterial species per niche differed significantly by queen age. Biological aging in queens was correlated with microbiota composition highlighting the relationship of microbiota with oxidative stress. Queens and workers shared many major gut bacterial species, but differ markedly in community structure and age succession. In stark contrast to aging workers, carbonyl accumulation in queens was significantly associated with increased Lactobacillus and Bifidobacterium and depletion of various Proteobacteria. Conclusions We present a model system linking changes in gut microbiota to diet and longevity, two of the most confounding variables in human microbiota research. As described for other model systems, metabolic changes associated with diet and host longevity correspond to the changing microbiota. The pattern of age-associated succession in the queen microbiota is largely the reverse of that demonstrated for workers. The guts of short-lived worker phenotypes are progressively dominated by three major Proteobacteria, but these same species were sparse or significantly depleted in long-lived queen phenotypes. More broadly, our results suggest that lifespan evolution formed the context for host-microbial interactions and age-related succession of honey bee microbiota.

decreases in core Bifidobacterium and Lactobacilllus, the same general results found in many 109 other microbiota models including insects and mammals [12,30]. Unlike workers, the queen does 110 not show greater antioxidant expression with age suggesting that antioxidant function is 111 performed differently or managed by her diet [9]. Vitellogenin hemolymph concentration, 112 constant royal jelly ingestion, and perhaps the microbiota contribute to antioxidant function in 113 long-lived queens.

115
While research on the worker microbiota has progressed rapidly, little is known of queens. Based 116 on a small sample size and 16S rRNA gene sequencing (amplicons) from whole guts, the queen 117 and worker microbiota differ in taxonomic membership and community structure [19,20,31]. 118 Unlike workers, the early queen gut seems dominated by two distinct species of 119 Acetobacteraceae, P. apium and an unnamed species referred to as "Alpha 2.1" [12]. Alpha 2.1 is 120 prevalent in guts of older workers, but P. apium occupies a variety of nutrition rich niches 121 associated with honey bees and thrives in the presence of royal jelly [7,32,33]. Capable of gut 122 colonization, P. apium is correlated with disease agents in adult bumblebees [34], and in honey 123 bees, implicated in poor worker health, increased mortality, worker gut dysbiosis, and strain-124 dependent effects on larval and pupal survival [12,29,33]. Often occuring with P. apium, 125 Lactobacillus kunkeei is prevalent/abundant in queens, but like P. apium, is also associated with 126 worker disease and dysbiosis [12,35]. L. kunkeei is also considered a hive (not a gut) bacterium 127 due to its association with fructose rich niches like honey and honey rich pollen storage [32,36]. 128 Thus, investigations of honey bee microbiota require a careful consideration of social and 129 functional context including host longevity, caste specificity, developmental stage, potential 130 refugia and transmission from nutrition related niches [22,23,37]. 131 Here we test the hypothesis that lifespan differences in a social insect model are associated with 132 age-based microbial succession. We sample known age queens from different backgrounds, and 133 compare our findings to the extensive preexisting characterization of known age workers. We 134 define the aging queen microbiota by deep sequencing four alimentary tract niches that differ in 135 many ways including physiological function, pH and oxygen exposure. To accompany each 136 amplicon library, we determine the absolute numbers of bacteria with qPCR. Finally, we 137 quantify protein oxidation in the fat body tissue of each queen to test the hypothesis that 138 biological age differs from chronological age, and that the accrual of oxidation products in aging 139 queens is associated with species -specific differences in the microbiota.

142
Queen sampling 143 Our sampling design distinguished environmental exposure from chronological age. We We sampled a total of 63 queens. Half (n= 32) of these queens were sourced from a large 151 migratory beekeeping operation based in southern California. Referred to as the "CA" source, The other half of our sampled queens, referred to as "AZ" queens (n = 31) were sourced from 161 two very different environmental and genetic backgrounds. We sampled "young" AZ queens (n 162 = 15) from the Carl Hayden Bee Research Center in Tucson Arizona. Delivered and installed 163 with 3000 young worker bees (package bees) in early May 2016, these Italian queens (Apis 164 mellifera ligustica) were exposed to varied pollen and nectar sources typical of the Sonoran 165 desert, but not intensive agriculture or bulk transportation events. Young AZ queens were 166 sampled in early October, 2016 at 5.7 months of age. In contrast, old AZ queens originated from 167 a Northern migratory beekeeping operation that raises Apis mellifera carnica queens. These All 63 queens were collected into sterile 2.0ml tubes and immediately frozen on dry ice and 177 stored at -20°C for DNA extraction. Queens were dissected under sterile conditions. Four tissue 178 types were extracted to typify the queen microbiota: mouth parts, midgut, ileum and rectum. were joined using the make.contigs command. After the reads were joined the first and last five 230 nucleotides were removed using the SED command in UNIX. Using the screen.seqs command 231 sequences were screened to remove ambiguous bases. Unique sequences were generated using 232 the unique.seqs command. A count file containing group information was generated using the 233 count.seqs command. Sequences were aligned to Silva SSUREF database (v102) using the 234 align.seqs command. Sequences not overlapping in the same region and columns not containing 235 data were removed using the filter.seqs command. Sequences were preclustered using the 236 pre.culster command. Chimeras were removed using UCHIME [42] and any sequences that were 237 not of known bacterial origin were removed using the remove.seqs command. All remaining 238 sequences were classified using the classify.seqs command. All unique sequences with one or 239 two members (single/doubletons) were removed using the AWK command in UNIX. A distance  honeybee-associated gut microbiota on NCBI (accessed July 2013). OTUs were generated using 245 the cluster command. Representative sequences for each OTU were generated using the 246 get.oturep command. To further confirm taxonomy, resulting representative sequences were 247 subject to a BLAST query using the NCBI nucleotide database. Diversity indices were 248 generated using the rarefaction.single and summary.single (alpha diversity) commands.

250
Statistical analysis 251 To examine the effect of community size we multiplied the proportional abundance of OTUs 252 returned by amplicon pyrosequencing by the total bacterial 16S rRNA gene copies determined 253 with qPCR for each individual queen and niche. All core bacterial genomes contain four 16S 254 rRNA gene copies except L. kunkeei (5), Bifidobacterium (2) and P. apium (1). Acetobacteraceae 255 Alpha 2.1 (copy number unknown) was designated a value of one, consistent with the copy 256 number of its closest relative, P. apium [44]. OTUs representing non-core diversity were 257 summed (Σ OTUs 10-500), corrected for community size and mean 16S rRNA gene copy  . We compared microbial community structure by chronological age and 265 source using a two-way factorial MANOVA and a post-hoc test to compare specific bacteria 266 across conditions. We compared absolute abundance of each bacterial taxon by age without 267 reference to source variation using the Wilcoxon rank sum test. Finally we examined the 268 relationship between carbonyl accumulation and the microbiota in various ways: 1) Using 269 DistLM we test whether the microbiota from each of four distinct tissues is significantly 270 associated with carbonyl accumulation in queens, 2) We examine carbonyl accumulation as a 271 covariate in three separate MANCOVA models, a two-way examining source and age, a one-272 way examining source, and a one-way examining age, 3) We calculate independent Pearson's 273 correlations between species-specific CLR scores and log transformed carbonyl data, and 4) We  and prevalence across many worker studies, and was combined as an "OTHER" category to 289 represent low abundance bacteria or signs of dysbiosis. As stated above for queens, we CLR 290 transformed relative abundance measures of workers, and performed a one-way MANOVA on 291 age to compare forager vs. in-hive bee gut microbiotas, calculating post-hoc differences between 292 specific bacterial groups. To compare the queen and worker results we transform our tissue-293 specific queen data to reflect relative abundance values predicted for the whole gut. Tissue 294 specific bacterial cell counts were used to normalize the relative occurrence of bacterial species 295 by queen tissue, then additively produce a single value that represents the expected result of 296 sequencing whole queen guts. These whole gut values highlight differences in abundance and 297 prevalence between workers and queens.

300
Next generation sequencing and qPCR 301 Next generation sequencing returned 7.2 million quality trimmed reads (400 bp assembled) 302 across the 252 libraries (63 queens X 4 niches). Read coverage was sufficient for all downstream 303 characterization and statistics (Table S1). The queen rectum was represented by 2.4 million reads 304 averaging 38K per library, the ileum by 1.9 million reads averaging 30K per library, the midgut  (Table S2). Species dominance in the queen increases with community size in the mouth, 322 midgut and ileum (Table S2). Extrapolating qPCR results to estimate absolute abundance, P. 323 apium, and L. kunkeei decrease in relative abundance approaching the rectum, while L. firm5, L. 324 firm4, Bifidobacterium and Alpha 2.1 increase (Table S3). S. alvi and G. apicola occur 325 sporadically at low (< 1%) relative abundance throughout all queen niches.

328
The two-way MANOVA performed for each of the four queen niches revealed significant 329 variation due to chronological age, source and interaction ( Table 1, Table S4). In the mouth, P.  (Table S3). In the rectum, 335 where L. firm5 represents the majority of total gut bacteria, Bifidobacterium abundance increased 336 with age while L. firm5 and both core Acetobacteraceae (P. apium and Alpha 2.1) decreased 337 (Fig. 2). Wilcoxon rank sum tests revealed significant differences by chronological age, many of 338 which agree with age-specific differences detected in the two-way MANOVA ( Table 1, Table   339 S4). proportional shifts among core gut bacteria (Fig. 1). All three core Firmicutes (Bifidobacterium, 345 L. firm5 and L. firm4) decreased significantly with age while "OTHER" bacteria, 346 Acetobacteraceae Alpha 2.1, and Bartonella apis increased significantly (Table S5). Of note, S. 347 alvi decreased in relative abundance with age but was borderline insignificant (p = 0.06). With 348 reference to these results, and other worker data sets in the literature we define four worker-349 specific gut species, all Proteobacteria, two queen-specific species, and four species shared by 350 longevity phenotypes (Fig. 3). In general, the species shared by longevity phenotypes are 351 particular to the rectum while the ileum species show fidelity by longevity phenotype.

353
Molecular age and the queen microbiota 354 We measured carbonyl accumulation in the queen fat body as a proxy for queen molecular age.

355
While some variation in carbonyl accumulation is due to genetics and background, difficult to 356 excrete waste products accumulate in a clock-like fashion with age. We found that chronological 357 age was strongly associated with carbonyl content in the fat body of the queen (Fig. 4). Carbonyl 358 accumulation differed by both age and source (Table S6). Examining all pairwise combinations, 359 only first year queens (CA1 and AZ1) did not differ in average carbonyl accumulation. In both 360 sets of young and old queens, chronological age did not agree with molecular age. In both age 361 classes, queens from the Imperial Valley of California (source CA) were chronologically 362 younger, but biologically older with greater carbonyl accumulation (Fig. 4).

364
To further explore the relationship of carbonyl accumulation with queen microbiota, we  (Table S8). Most notably, throughout the 377 gut Bifidobacterium is correlated significantly with the accumulation of carbonyl in abdominal 378 fat body tissue (Fig. 5). Although at similar abundance in chronologically old and young queens, 379 L. firm5 abundance was also correlated strongly and positively with carbonyl accumulation.

380
Although rare throughout the queen gut, an undescribed Burkholdariales; Delftia showed the 381 strongest negative relationship with carbonyl content, decreasing dramatically with age, and 382 varying by source (Table 1).

384
To better visualize variation associated with biological age in the queen microbiota, we 385 performed PCA analysis using centered log ratios from the top 9 OTUs and associated carbonyl 386 values from the fat body of each queen (Fig. 5, Table S9). Across each niche the first two 387 principle components explained approximately 50% of the variation in log ratio abundance 388 scores. Because the queen microbiota has shallow, deep and noisy structure, the third and fourth 389 principle components for each niche explained an average of 14% and 8% respectively (Table   390   S8). Although only 50% of the variation is presented in the two dimensional PCAs, a strong and 391 consistent separation of two queen cohorts is realized in every niche; young Arizona (AZ1) and 392 old California (CA2). In each niche, the carbonyl vector indicates CA2 as the oldest, and AZ1 as 393 the youngest cohort, consistent with determinations of molecular (biological) age.

395
We examined microbiota correlations using SparCC, an approach that incorporates the structure  (Table S10). We note that SparCC results are unreliable when OTU sparsity exceeds 70% 400 zero values but robust to communities with a low effective number of species (see Table 1). The Bifidobacterium. The strongest negative correlation occurs between L. firm5 and P. apium, the 405 two most abundant ileum species (Table S10). With more detailed investigation, age-specific 406 Pearson's correlations on log transformed absolute abundance shows that as queens age, the

412
We show that host phenotypes with extreme longevity differences support gut microbiotas that 413 age differently (Fig. 1). Because long and short-lived phenotypes are produced from the same 414 genotype, microbiota establishment and age-associated changes likely reflect host gene Longevity phenotypes differ in core membership 430 Microbiotas of long and short-lived phenotypes differ markedly in core bacterial membership 431 sharing four of ten species, with six species showing strong phenotype-specificity (Fig. 3). The  The queen microbiota improves with age 464 Our results are consistent with the body of work detailing molecular aging and oxidative stress 465 in queens, workers and social insects in general [8-10,65]. Results from the queen carbonyl 466 assay demonstrate that queens accrue oxidative damage with age (Fig. 4), and that chronological 467 age can differ significantly from biological age possibly due to environmental differences 468 including climate, nutrition, toxins and other landscape variables. Despite similar signs of 469 biological aging in both queens and workers, the gut microbiota of older queens seems to reflect 470 a refined structure with greater efficiency. It's unlikely that queens ever develop a senescence 471 physiology and associated microbiota as seen in workers. Under natural conditions, queens 472 accrue molecular damage associated with aging, but are not allowed to grow "old" because 473 fecundity is critical to colony survival, and workers routinely replace substandard queens [66]. 474 With increased oxidative damage, gram positive bacteria decrease in workers [19] but increase 475 in queens (Fig 1). Of note, core Lactobacillus and Bifidobacterium in queens show greater 476 correspondence with biological than chronological age suggesting that these species may track or 477 signal host physiology (Table S7). Consistent with decreased antioxidant expression and less 478 ROS generation in queens [9,65] the bacteria that increase with queen age do not rely on oxygen, 479 but generate continuous fermentative metabolism in the queen hindgut (Fig. 2). In turn, these . We found that bacterial communities implicated in butyrate production were 489 diminished in aging workers but seemingly enhanced in aging queens. Better explained by 490 biological than chronological age (Fig. 5), Bifidobacterium increases significantly with age in the 491 queen midgut, ileum and rectum (Table S7). Moreover, L. kunkeei increases in the midgut and 492 ileum, while community changes in the ileum favor the persistence of L. firm5 (Table 1)  We compared the gut microbiota of young in-hive bees to older foragers within and among 502 studies. Foraging is the last functional role workers serve before death. But as a group, both in-503 hive bees and foragers can range greatly in chronological age and environmental exposure [69]. 504 Also, comparing across next generation sequencing studies can be misleading due to differences 505 in methodology like primer choice or analysis pipeline [70]. Despite these and other sources of 506 potential error, we found that the worker gut microbiota ages in a highly predictable fashion, 507 becoming depleted of core hindgut firmicutes including Bifidobacterium (Fig.1). Of seven 508 available forager studies, three used the same methods to sequence both foragers and nurses 509 [19,31,63], and we used these studies as a point of reference for examining worker aging. We 510 analyzed the largest and most variable of these three data sets [19] defining six significant 511 differences in microbiota between young and old workers (Fig. 1). The collective results from six 512 of seven studies are largely in agreement and suggest that age-associated shifts in worker 513 microbiota are strongly predictable at the level of species despite study particulars [19,31,60-514 64]. The changes we report in figure 1 [19] represent a functional change from a fermentative to 515 proteolytic hindgut environment, involving significant shifts in core bacterial structure. Alpha 516 2.1 increases in all 7 studies, B. apis and "Other" bacteria in 6 of 7. One to three major core 517 hindgut Firmicutes are depleted significantly in 6 of 7 studies, while studies were more variable 518 concerning shifts of S. alvi and G. apicola, the species pair that dominates the worker ileum.

520
Early gut succession 521 Similar to workers [37], the rectum of the mature queen contains 84% of the total bacteria found 522 in the queen gut (Fig. 2). On average, a whole gut sample from a mature laying queen would be 523 highly biased toward rectum species, dominated by Lactobacillus firm5 (Fig. 3). In contrast, Queen niche breadth 538 In queens, the occurrence patterns and numerical dominance of P. apium in the mouth and 539 midgut, and L. firm5 in the ileum and rectum suggests that the extended gut structure is 540 important for host function (Fig. 3). The taxonomic shift at the pylorus demarcates a steep products, supporting hindgut bacterial strains distinct from those found in workers.

551
It is mostly unknown why queens can resist many worker diseases and vice-versa. Early queen 552 death has become more common [66,72], and defining disease states in queens will rely in part 553 on the structure and function of native gut bacteria [12]. Although rare throughout the gut, the apium on their mouths and L. kunkeei in their midguts that may accrue with age and/or improve 568 queen hygiene promoting queen survival (Fig. 2). Pollen exposure and consumption may render  The evolution of host behavior to mechanically concentrate nectar sugars via evaporation was 581 likely a key innovation producing strong selection for these two osmotolerant symbionts. . In contrast, the mature queen ileum is dominated by Lactobacillus firm5 that co-586 occurs with lesser amounts of core gut bacteria P. apium and L. kunkeei (Fig 1). Reciprocally, 587 worker ileum bacteria S. alvi and G. apicola are found at similarly low levels in the queen ileum 588 and show sporadic abundance in the queen. These symmetrical occurrence patterns suggest 589 antagonistic co-evolution of caste-specific gut bacteria, a hypothesis consistent with host age 590 phenotype and development-specific pathogen strategies. Patterns of species co-occurrence suggest selection pressure for honey bee gut bacteria to co-607 exist with other bacteria in a biofilm encouraging competition and co-evolution (Fig. 4). This shifts towards the fermentative metabolism of well-known gram positive species, while the more 625 rapidly aging worker is progressively depleted of these same species. Given the spectrum of 626 influence of gut microbiota on worker physiology, we suggest that the queen microbiota serves a 627 similarly critical role in host signaling and protection. Separate evolutionary trajectories for 628 caste-specific gut bacteria reflect overt differences in diet and longevity between workers and 629 queens. This trajectory appears to have tracked division of labor evolution, perhaps involving 630 key innovations like nectar concentration to produce honey, and the production of royal jelly in 631 worker hypopharyngeal glands. Once considered bacteria associated with worker gut dysbiosis 632 and larval nutrition, L. kunkeei and P. apium must now be understood as core gut bacteria of Apis 633 mellifera queens. Our results suggest that these two species occupy a functional niche in the 634 queen mouth, midgut and ileum. The co-occurrence and correlational abundance of multiple core 635 species throughout the honey bee system suggest syntrophic relationships are commonplace.

636
More generally, our study highlights the importance of controlled temporal and tissue-specific 637 data to understand the total diversity and function of the honey bee microbiome.

859
C Average percent change in bacterial cell number with age. We note that cell number loss cannot exceed 100%.

860
D Comparing species-specific bacterial cell number by chronological age only (Table S4).

861
E Independent variables are queen age and background (df = 3, 59). Reports only F-values for chronological age effects

863
* Significant interaction effect of age and background detected by the 2-way MANOVA (Table S8).  Table S3 for 880 absolute abundance). The 4x4 panel displays the top 9 most abundant OTUs by niche, age and 881 source. Black represents "diversity abundance", the summation of OTUs 10-500. Old queens in 882 the upper two rows are 16-18 months of age and young queens in the bottom two rows are aged 883 4.5-5.7 months (Fig. 4).   Workers are whole gut samples from Kwong et al. [19]. Queen data was normalized by tissue-894 specific community size to reflect relative abundance values expected from sampling whole guts.

895
The red bars represent average abundance, black bars are prevalence defined at ≥ 0.5% relative 896 abundance, and the bar apex is prevalence defined as 2 or more reads per gut library. We did not 897 detect F. perrara* in any of the four sampled queen alimentary tract niches.   (Table S9).