3 RESULTS

3.1 16S rRNA gene sequencing

Next-generation sequencing returned 11.5 million quality trimmed reads (455 bp assembled) across 292 libraries. The worker midguts represented about 4.4 million reads averaging 31 K per library. The worker ileums represented about 7.1 million read averaging 49 K per library. A total of 1230 OTUs were resolved at 97% similarity. The top 10 OTUs and a 11th group consisting of ‘Other’ (Σ OTUs 11-1230), represented 95% and 5% of the total reads, respectively. The most abundant OTUs according to raw amplicon read totals were Gilliamella apicola (23.1%),Snodgrassella alvi (17.6%), Frischella perrara (13.2%),Apilactobacillus apis (Firm5-1, 12.7%),Apilactobacillus kunkeei (7.3%), Apilactobacillus Firm4 (7.1%), Apilactobacillus Firm5 (Firm5-2, 7.1%),Bifidobacterium sp. (3.8%), Fructobacillus fructosus(2.5%), and an Enterobacter sp. (0.8%) (Figure 2A). This is largely consistent with previously described worker microbiota analysis (Kwong & Moran, 2016). Our analyses clustered ApilactobacillusFirm-5 as two distinct clades, Apilactobacillus apis(“Firm5-1”) which favors the ileum and a second (“Firm5-2”) cluster that was more abundant the midgut (Figure 2A). Firm5-2 BLAST hits gave sequences identities to Lactobacillus panisapium, Lactobacillus kimbladii, and Lactobacillus helsingborgensis .
Our MC and SCC age-right nurses gut microbiota were compared and did not differ statistically (Table S2). Therefore, we added MC nurses to 7 and 13-day old SCC nurses for all downstream statistics.

3.2 Microbiota associations with age and task

The two-way MANOVA performed on midguts and ileums revealed significant variation due to age and behavioral task, but not as an interaction factor (Table 1). In the midgut, the MANOVA revealed significance for eight of the OTUs analyzed. The independent variables (IV) of age (F value 3.72, Pr > F = 0.0001) and behavioral task (F Value 5.56, Pr > F = 0.0001) were significant for six OTU each. Specifically, G. apicola , was more abundant in age-right nurses versus their same age PFs. G. apicola also increased with age; in 27-day old age-right foragers and 19-day overage nurses had the largest proportions and the largest microbiotas based on 16S rRNA gene copies (Figure 3A). S. alvi and Enterobacter sp. increased with age regardless of behavioral task, while A. kunkeei and F. fructosus showed larger relative abundances in nurses rather than foragers and were strongly correlated (Table S3: Spearman’s ρ, rs= 0.90, p = 0.0002). Bifidobacterium sp. , Firm4, and Firm5-2 were highly correlated across all samples (Table S3: Spearman’s ρ, rs= 0.90, p = 0.0002), but had higher relative abundance in PFs versus age-right nurses as well as increasing with age. PCA on CLR adjusted midgut microbiota shows tight grouping for behavioral task, with 31.6% explanation for the first component and 15.1% for the 2nd component. Grouping is also consistent on microbial based age-association (Figure S1). Patterns predicted by MANOVA are represented in the PCA withA. kunkeei and F. fructosus being more associated with nurses and the highly correlated OTUs, Bifidobacterium sp , Firm4, and Firm5-2 found in the guts of foragers.
A visual inspection of the ileum’s relative abundances shows remarkable stability across age and behavioral task compared to the midgut (Figure 2A). The size of the ileum microbiota was mostly consistent across treatment groups (Figure 3A). Nevertheless, the two-way MANOVA revealed significant variation due to age (F value 3.17, Pr > F = 0.0001) and behavioral task (F value 4.74, Pr > F = 0.0001). There were five OTU significant in the model, five of which were best explained by age and four by behavioral task. The bacteriaG. apicola, S. alvi, Firm4, and Firm5-2 were explained by both factors in the model. Firm4 and Firm5-2 were also very strongly correlated across all samples (Table S4: Spearman’s ρ, rs = 0.82, p = 0.0006) A PCA on CLR adjusted ileum microbiota shows overlap for behavioral task, with 23.7% for the first component and 19.5% for the 2nd component. Age-associated grouping is also consistent with the MANOVA with less overt differences in the ileum relative to the midgut (Figure S2).
Wilcoxon on absolute abundance were largely in congruence with the significant results of the two-way MANOVA (Table S5 and S6). There was also uniformity with gut microbiota size; 27-day age-right foragers had significantly more cell counts for many of the top OTUs in the midgut (Figure 3A). 16S rRNA gene copies increased significantly with age in the midgut but the ileum microbiota remained the same size over time. The fungal loads were more similar between sample types and tissues but bifurcated by age and task in the midguts; the oldest nurses and 27-day foragers had the greatest fungal loads (Figure 3B). The alpha diversity metric for richness was significant for age and sample type (p -value < 0.0001), while beta diversity was significant for task only ([PERMANOVA] F-value: 5.3629; R-squared: 0.01991; p -value < 0.02) (Figure S5).
3.3 Gene expression associations
Two-way ANOVA found an interaction effect between age and behavioral task for all immune gene expression in the fat body, except for the AMPs apidaecin and hymenoptaecin (Table S7). Nonetheless, both AMPs were significant for age or behavioral task independently.
Two PCAs reveal patterns in immune gene expression by behavioral task and/or age (Figure 4). The first principal component (29.1%) maximizes the variance of immune gene expression, with gene vectors along the horizontal explained more by behavioral task (Figure 4A). In contrast, genes that are parallel to the second component (20.1%) are more closely associated with age (Figure 4B). Significant differences between sample types dichotomize by behavioral task. Vitellogenin is highly expressed in all age nurses and downregulated in age-matched and age-right foragers (Figure 5A). Hymenoptaecin and catalase had increased expression in foragers (Figure 5B-C). Apidaecin, CuZnSOD, and GST-1 followed age-associated gene expression patterns with age-right 27-day foragers having the highest transcript levels (Figure 5F-H). DSCAM and MnSOD had expression patterns that were affected by task and age (Figure 5D-E). DSCAM was expressed higher in age-matched PFs versus nurses but decreased in both with age. MnSOD was highest in age-right nurses compared to age-matched foragers, and decreased as the PF aged. MRJP2 was depressed in older nurses but also appears to increase with age while PPO was expressed at the same levels across all sample types (Figure S6).
Midgut tissue showed the most variation in microbial composition (Figure 2A) and size (Figure 3A). Therefore, we looked at correlations between the midgut bacterial absolute abundance with fat body gene expression to investigate potential relationships between the aging midgut and overall immune health (Table S8). In PFs only, we found a positive correlation of CuZnSOD with F. fructosus (rs = 0.42,p < 0.008). There were no significant correlations in nurses, but when considering all samples, including age-right 27-day foragers, we found many positive relationships with CuZnSOD; G. apicola (rs = 0.32, p < 0.0001), Firm5-1 (rs = 0.29, p < 0.002),F. fructosus (rs = 0.30, p < 0.002), A. kunkeei (rs = 0.29, p< 0.002), S. alvi (rs = 0.27, p< 0.005), and Bifidobacterium (rs = 0.24, p < 0.01). Hymenoptaecin gene expression correlated positively with Bifidobacterium (rs = 0.2, p < 0.006) and Firm5-2 (rs = 0.23,p < 0.02). There were two negative correlations found with all samples, MnSOD with Firm5-2 (rs = –0.23,p < 0.02) and DSCAM with F. fructosus(rs = –0.32, p < 0.05).

3.4 Carbonyl contents of abdomen

Oxidative damage in the abdomen was measured by carbonyl assay of the fat body and attached abdominal sclerites. Foragers had accumulated more oxidized proteins than age-matched nurses, a finding that became statistically meaningful after 19-days (Figure 6). 19-day foragers and nurses had also accrued significantly more carbonylation than age-right 27-day foragers. As predicted Carbonyl contents decreased with increased Vg expression (rs = -0.28 p < 0.0006). To consider whether fat body gene expression or midgut bacterial absolute abundance had a relationship with carbonyl content accumulation we ran a Spearman’s correlation analysis (Table S8). Considering all samples, carbonyl contents were negatively associated with the abundance ofG. apicola (rs = –0.26, p < 0.0071). The abundance of S. alvi and G. apicola was positively correlated in workers with natural ontogeny (nurses; rs = 0.55, p < 0.01), but showed no relationship in precocious foragers.