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.