2 METHODS

2.1 Colony Manipulations and Sampling

In June 2019 at the USDA-ARS Carl Hayden Bee Research Center in Tucson Arizona, closed brood frames were sourced from 30 honey bee (Apis meliffera, linguistica) hives. Frames were incubated overnight (30°C, 75% relative humidity, and 24h dark cycle) and newly emerged adult workers were collected and combined into a mixed cohort to implement into experiment. We utilized a “single-cohort colony” (SCC) design, which was previously shown to uncouple behavioral task from chronological age (Robinson et al., 1989) (Figure 1A). To construct the two SCCs, newly emerged workers (4500 and 3500) were added to a small hive box, each containing a naturally mated queen, one frame with pollen and honey, and one frame with eggs and open brood. Bees assigned to each SCC were less than 24hours old and differentiated into separate behavioral tasks, i.e., nurses and atypically, precocious (young) foragers (PFs). Additionally, a marked cohort (MC) with 4100 newly emerged workers was constructed to serve as a control for sampling normal ontogeny (Figure 1A). Newly emerged workers were marked with paint on their thorax and transferred into a healthy double-deep colony free from visible signs of disease. On day 6 of the experiment, SCC PFs we observed with corbicular pollen loads were marked with paint on their thorax.
At the peak of foraging activity in the summer, a workers lifespan is approximately 30 days (Fluri et al., 1982). Over the course of the experiment, we sampled foragers and nurses from the SCCs at 7, 13, and 19-days old (Figure 1B). By design, SCC nurses and foragers performed the same behavioral tasks for the duration of the experiment. Nurses were identified by observing the brood nest and sampling bees that spent 3 seconds with their head in a cell containing brood. Part way through the experiment, we replaced the brood frame in the SCC’s to ensure no newly emerged workers replace the current nurses. For the MC sampling, we sampled age-right nurses at 7 and 13-days old and age-right foragers at 27-days based on well-established honey bee worker ontogeny (Seeley, 1982). Thus all samples could be categorized as normal (age-right) or atypical (PFs and overage nurses) ontogeny (Figure 1C). Sampled bees were collected with sterile soft forceps, snap frozen with dry ice, and stored in -80°C for processing.
Worker dissections occurred under sterile conditions. The sting was discarded and the fore and hindguts were removed from the abdomen. Gut tissues were dissected in 70% EtOH to wash and separate the midgut and ileum before being added to a bead-beating tubes with 0.2 g of 0.1-mm silica beads and 600 µL of 1X TE buffer. The abdominal fat body and attached dorsal sclerites were retained in similar bead beating tubes as a single unit for gene expression and protein oxidation (carbonyl) assay.

2.2 DNA extraction and qPCR

In preparation for DNA/RNA extractions, samples were bead-beaten for 2 min at 30-s intervals and centrifuged to recover the supernatant. Midgut and ileum DNA were extracted according to Thermo Scientific GeneJET Genomic DNA Purification Kit manufacturer’s instructions. Fat body supernatant was split 300 µL for the carbonyl assay and 300 µL for RNA extraction with Thermo Scientific GeneJET RNA Purification Kit according to manufacturer’s instructions. The extracted fat body RNA was converted into cDNA with Thermo Scientific RevertAid First Strand cDNA Synthesis Kit according to manufacturer’s instructions.
We quantified total bacterial abundance for the midgut and ileum with a real-time PCR (qPCR) assay of bacterial 16S and fungal 18S rRNA gene copies (C. M. Liu, Aziz, et al., 2012; C. M. Liu, Kachur, et al., 2012). This assay provides significantly broader coverage than previously reported universal quantification assays. The bacterial 16S gene template was amplified using forward primer 27F (5’-AGAGTTTGATCCCTCAG -3’) and reverse primer 1522R (5’- AAGGAGGTGATCCAGCCGCA -3’). The fungal 18S gene template was amplified using forward primer PanFungal_18S_F (5’- GGRAAACTCACCAGGTCCAG -3’) and reverse primer PanFungal_18S_R (5’-GSWCTATCCCCAKCACGA-3’). Invitrogen’s pCR™2.1 TOPO™ was used to produce plasmid vectors, which were then transformed into DH5α™ cells. Successfully transformed colonies were selected and grown overnight in broth. Plasmid DNA was purified using the Thermo Scientific GeneJET Plasmid Miniprep Kit according to manufacturer’s instructions. The purified plasmid cells were measured using an Implen nanophotometer P300 to calculate 10-fold serial dilutions as the standards for qPCR quantification.
Quantitative PCRs for 16S and 18S were carried out in triplicate on a BioRad CFX96 thermocycler in 12 µL reactions containing 6 µL of New England Biolabs - Luna® Universal Probe qPCR Master Mix, 0.5 µL forward primer, 0.5 µL reverse primer, 3 µL of H2O and 2 µL of DNA template. The cycling conditions were 95 °C for 3 min followed by 40 cycles of 95 °C for 10 s and 60 °C for 60 s. The qPCR results were expressed as the total number of 16S and 18S rRNA gene copies per DNA extraction (100 µL volume elution).

2.3 Immune Gene Expression

The fat body is a main metabolic tissue of the honey bee and is functionally analogous to vertebrate liver or adipose tissue (Y. Liu et al., 2009). Comparisons of fat body gene expression can relay information on immunocompetence and overall health. Fat body immune and oxidative stress gene expression was quantified for vitellogenin(vg ), apidaecin , catalase , copper-zinc superoxide dismutase (CuZnSOD ), down syndrome cell adhesion molecule (DSCAM ), glutathione S-transferase-1 (GST-1 ),hymenoptaecin , manganese superoxide dismutase (MnSOD ), major royal jelly protein 2 (mrjp2 ), and prophenoloxidase (PPO ).
Quantitative PCR reactions for immune gene expression were performed in triplicate as follows: initial denaturation at 95 °C for 5 minutes; 40 cycles with denaturation at 95 °C for 15s; and a primer-pair-specific annealing and extension temperature for 30 seconds (Table S1). To confirm the absence of contaminating genomic DNA and primer dimers in the qPCR assay, we monitored amplification and melting curves. Relative gene expression was determined based on standardized Ct values (Δ Ct) (Livak & Schmittgen, 2001) using actin as a reference gene. Actin is constitutively expressed in different honey bee tissues and has been previously established as an effective control for measuring gene expression in adult workers (Evans et al., 2013).

2.4 Carbonyl Assay

To measure protein damaged by oxidative stress, we quantified the accumulation of protein carbonyl groups (Reznick & Packer, 1994). The method for measuring carbonyl content in honey bees has been previously established (Anderson et al., 2018; Hsieh & Hsu, 2011; Williams, Roberts, & Elekonich, 2008). To determine carbonyl content of fat body homogenates, we used Protein Carbonyl Content Assay Kit (MAK094; Sigma-Aldrich). Briefly, samples were treated with a 10 mg/ml streptozocin solution and incubated for 15min to precipitate nucleic acids. Keeping the supernatant, 2,4-dinitrophenylhydrazones (DNPH) was added to samples to form stable dinitrophnyl hydrozone adducts. Derivatized proteins were precipitated with trichloroacetic acid and were followed by three successive ice-cold acetone washes. Samples were resuspended in 100 μl of 6 M guanidine (pH 2.3). The total protein concentration of each sample was measured using a PierceTM BCA Protein Assay Kit (Smith et al., 1985). Protein oxidation was expressed as nanomoles of carbonyl groups per mg of protein.

2.5 16S rRNA gene community analysis

16S rRNA gene sequences were processed using MOTHUR v.1.44.3 (Schloss et al., 2009). Forward and reverse reads were joined using the make.contigs command. After the reads were joined, the first and last five nucleotides were removed using the SED command in UNIX. Sequences were screened to remove ambiguous bases, using the screen.seqs command. Unique sequences were generated using the unique.seqs command. A count file containing group information was generated using the count.seqs command. Sequences were aligned to Silva SSUREF database (v102) using the align.seqs command. Sequences were filtered to remove overhands at both ends and gaps using filter.seqs. The unique.seqs command was ran again to remove new redundancies from filtering. A precluster step using pre.cluster was performed. Chimeras were removed using chimera.uchime command (Edgar, Haas, Clemente, Quince, & Knight, 2011). Sequences were classified with a bayesian classifier using the Silva database using classify.seqs command. Sequences that were not bacterial origin were removed using the remove.seqs command. All unique sequences with one or two members (single/doubletons) were removed using the AWK command in UNIX. A distance matrix was constructed for the aligned sequences using the dist.seqs command. Sequences were classified at the unique level with the Ribosomal Database Project (RDP) Naive Bayesian Classifier (Wang, Garrity, Tiedje, & Cole, 2007) using a manually constructed training set containing sequences sourced from the greengenes 16S rRNA database (version gg_13_5_99 accessed May 2013), the RDP version 9 training set, and all full-length honeybee-associated gut microbiota on NCBI (accessed July 2013). Operational taxonomic units (OTUs) were generated at the 97% species using the cluster command. Representative sequences for each OTU were generated using the get.oturep command. To further confirm taxonomy, resulting representative sequences were subject to a BLAST query using the NCBI nucleotide database.

2.6 Statistical Analysis

We evaluated our microbial datasets as both relative and absolute abundance for downstream analyses. To examine the effect of community size, the top ten OTUs and a sum of remaining OTUs were used for downstream analysis. It should be noted that OTU 1-10 accounted for 95% of all sequences, the 11th group consisted of ‘Other’ (Σ OTUs 11-1230) accounted for the remaining 5%. OTUs were normalized by 16S rRNA gene copy number via ribosomal RNA operons (rrn) database (Stoddard, Smith, Hein, Roller, & Schmidt, 2015) and total bacterial 16S rRNA gene copies from qPCR prior to analysis. In this case, qPCR-normalized abundance is extrapolated from relative abundance of amplicons, so remains compositional. To allow the use of parametric multivariate analyses (Pearson, 1897), we converted the qPCR-normalized bacterial abundances to ratios among all OTUs (Gloor & Reid, 2016) using the software CoDaPack’s centered log ratio (CLR) transformation (Comas & Thió-Henestrosa, 2011). After conversion, nearly all OTUs followed normal distributions. A two-way MANOVA was performed on log-converted relative abundances which allows for comparing taxa between dependent (OTU 1-11) and independent (age, behavioral task, and age*behavioral task) variables. As an additional measure, we used Pillai’s Trace test statistics which are robust in terms of violations of multivariate normality and homogeneity of covariance. To account for multiple comparisons, a False Discovery Rate (FDR) was conducted. We also performed principle component analysis (PCA) on CLR scores from OTUs 1–11, plotting the relationship of bacterial community composition with behavioral task and age-associated succession. To determine differences in absolute abundance of the microbial communities, we used Wilcoxon rank sum tests corrected for multiple comparison with Steel-Dwass. Absolute abundance was used to determine correlations between OTUs using Spearman’s ρ which was further corrected by FDR for multiple comparisons. MicrobiomeAnalyst web interface was used to analyze bacterial community diversity and composition (Chong, Liu, Zhou, & Xia, 2020). Alpha diversity metrics, Observed and Shannon were used to measure the richness and evenness using Mann-Whitney or Kruskal-Wallis statistical tests where appropriate. Beta diversity was used to compare the microbial communities with Bray-Curtis dissimilarity between each sample pair. The ordination was plot with non-metric multidimensional scaling (NMDS) using PERMANOVA as the statistical test for significance.
We evaluated bacterial and fungal copy numbers by sample type (age and task combined, i.e., 7-day nurse and 7-day forager) using one-way ANOVA with Tukey’s HSD post-hoc. Gene expression were log10 transformed to normalize variation and analyzed by sample type using TWO WAY one-way ANOVA with Tukey’s HSD post-hoc. PCAs of normalized gene expression was used to plot the relationship of immunity and oxidative stress with behavioral task and age. Carbonyl content significance was evaluated by sample type using TWO WAY one-way ANOVA with Tukey’s HSD post-hoc. In order to verify MC as a control following normal ontogeny we compared SCC age-right nurse microbiota and gene expression with the MC using Wilcoxon 2-sample t-test with FDR corrections. Multivariate analyses were conducted on OTUs with gene expression and carbonyl contents using Spearman’s ρ to find significant correlations after correcting for multiple comparisons with an FDR. All analyses were conducted in JMP_ v14.3.0 (JMP_ 1989–2007) and/or SAS_ v9.4 (Institute, 2015). We considered values of P < 0.05 statistically significant.