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.