Honey bee (Apis mellifera) queen quality: Host-microbial transcriptomes exploring the influence of age and hindgut symbiont Commensalibacter melissae

Understanding the biological mechanisms underlying extreme lifespan variation within species remains a fundamental challenge in aging research. Here, we investigated the role of gut microbiota and age in honey bee (Apis mellifera) queens combining metagenomics and transcriptomics. Analysis of 40 queen hindguts revealed that Commensalibacter melissae (Alpha 2.1) relative abundance was significantly higher in young queens compared to old queens. Using queens with the highest and lowest C. melissae relative abundance, RNA sequencing identified 1,451 differentially expressed genes associated with C. melissae abundance, twice the number associated with age alone (719 genes). Queens with high C. melissae abundance showed distinct transcriptional profiles related to stress response, protein homeostasis, and longevity-regulating pathways, particularly genes involved in oxidative stress response and cellular maintenance. Our analysis revealed complex relationships between age, C. melissae abundance, and gene expression patterns, suggesting that multiple interacting factors contribute to queen quality. These findings contribute to our understanding of host-microbe interactions in honey bee queens and highlight the intricate relationship between gut microbiota composition and host physiology in honey bees. 3.1 Queen Sampling Queens were sourced from two locations: the USDA-ARS Carl Hayden Bee Research Center in Arizona and a commercial beekeeping operation in Illinois. Arizona queens were all established and had survived more than one year in their colonies. Illinois queens represented a mixed-age population, consisting of both newly introduced queens (requeened in spring 2023, April-May) and older queens that had not been requeened that spring, though natural supersedure could not be ruled out. All queens were sampled in June 2023 from robust double-deep colonies. The selected colonies were highly productive with strong populations, often reaching space limitations for both brood production and honey storage. To ensure sample quality and consistency, we specifically selected queens from thriving colonies and deliberately excluded any colonies showing signs of queen failure or irregular egg-laying patterns. Forty queens in total were collected into sterile 2.0-ml tubes and immediately frozen on dry ice and stored at −80 °C for nucleic acid extraction. 3.2 Dissections, carbonyl protein oxidation assay, and nucleic acid extractions Queens were pinned through the thorax in 70% ethanol to wash and aid in dissection. Micro-dissection scissors were used to cut through the sides of the abdomen to access the digestive tract. The entire digestive tract was removed and floated in ethanol to manually separate the gut tissues with dissection tweezers. C. melissae is most abundant in the hindgut, so we targeted the ileum and rectum together for analysis. The abdominal fat body was extracted as a single unit for use in gene expression and protein oxidation assays to assess biological aging. Queen fat body, ileum, and rectum tissues were each bead-beaten separately in 1X TE buffer for 2 min at 30-second intervals and centrifuged at 30 seconds at 3000 rcf to recover the supernatant. The fat body supernatant fraction was used in a protein oxidation assay to quantify the accumulation of protein carbonyl groups associated with oxidative stress and aging (Reznick and Packer, 1994), as in Copeland et al., 2022b. Protein oxidation was expressed as nanomoles of carbonyl groups per mg of protein. To extract nucleic acids (DNA and RNA simultaneously) we used Qiagen AllPrep PowerViral DNA/RNA Kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol and methodology also reported in Copeland et al., 2022b. 20L of ileum and 20L rectum elution from the extractions were then pooled into hindgut samples for downstream sequencing and analysis. Using purified total DNA, we quantified total bacterial abundance for the hindgut using a quantitative PCR (qPCR) assay of the 16 rRNA gene (Liu et al., 2012). We created a standard curve using a 10-fold serial dilution series of a plasmid standard containing a full-length Escherichia coli 16S rRNA gene sequence. We amplified a 466 bp fragment in the V3–V4 region of the 16S rRNA gene using universal primer pair (5′-CCTACGGGDGGCWGCA-3′ and 5′-GGACTA CHVGGGTMTCTAATC-3′). PCR reactions were performed in triplicate on a BioRad CFX96 (Biorad, Hercules, California, US) as follows: 12 μl reactions containing 9 μl of iTaq Universal SYBR Green Supermix (BioRad, Hercules, California, US), 0.5 μl forward primer, 0.5 μl reverse primer, 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 rRNA gene copies per DNA extraction (200 μl volume elution). 3.3 Sequencing DNA from gut tissues was amplified in a single step procedure to amplify full-length 16S rRNA (V1-V9) using degenerate primers 27F (GCATC/barcode/AGRGTTYGATYMTGGCTCAG) and 1492R (GCATC/barcode/RGYTACCTTGTTACGACTT). PCR was performed with Q5 2× Hot Start High-Fidelity Master Mix (New England Biolabs) using the following conditions: 98℃ 30s; 98℃ 10s, 55℃ 30s, and 72℃ 2 min for 22 cycles; final extension of 72℃ for 10 min. Reactions (30 µL) were performed following manufacturer recommended master mix concentrations, with primer concentrations of 250 nM. Positive (ZymoBiomics Microbial Community DNA Standard; Zymo Research) and negative, non-template controls were included as process controls. After PCR amplification, target amplicons were purified from residual primers and primer-dimer using an AMPure bead cleanup and DNA concentrations were determined using a Qubit fluorometer (Thermo Fisher). Samples were then pooled (~ 3 ng per sample) and were prepared for sequencing by generating a SMRTbell library with a Pacific Biosciences SMRTbell prep kit 2.0 using manufacturer suggested inputs and procedures. Amplicons were sequenced on a single Pacific Biosciences 8M SMRT Cell on a PacBio Sequel IIe (Pacific Biosciences) at USDA-ARS PBARC (Hilo, HI). After sequencing, circular consensus sequences from the subreads were obtained using the SMRTLink v8.0 software. RNA from extracted samples was processed for sequencing by depleting the rRNA from samples using a RiboFree cDNA Kit (Zymo Research). Samples were pooled in equimolar concentrations and libraries were prepared for sequencing using an Adept Rapid PCR-Plus Kit (Element Biosciences). Libraries were sequenced on an Element Biosciences AVITI sequencer using an AVITI 2×150 Cloudbreak High Output sequencing kit at USDA-ARS PBARC. One queen sample dropped out of sequencing due to insufficient cDNA yield, leading to a failure in library preparation. Consequently, this sample was excluded from downstream RNA-seq analyses. All paired-end raw reads were filtered and trimmed using Trimmomatic v0.38 (Bolger et al., 2014) with the following parameters: LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:75. FASTQC was used to ensure quality control (Babraham Bioinformatics - FastQC A Quality Control tool for High Throughput Sequence Data). Kraken2 was used to map and split reads to the Apis mellifera genome assembly Amel_HAv3.1 (PRJNA471592) and the gut symbiont Commensalibacter melissae (PRJNA495947) (Wood et al., 2019). A. mellifera reads were mapped against A. mellifera genome using STAR v2.7.10b (Dobin et al., 2013) and C. melissae reads were mapped to C. melissae genome using Bowtie 2 v2.5.2 (Langmead and Salzberg, 2012). Gene counts were obtained using Subread v2.0.4 package FeatureCounts (Liao et al., 2013). Genes with read counts below 4 were removed, and genes with variance less than 15% across samples were filtered out. Final counts were normalized by employing a log2-counts per million (logCPM) transformation.