While honey bees play a vital role in global crop production, they face increasing exposure to xenobiotic chemicals during commercial pollination. Multidrug-resistance (MDR)-type ATP-binding cassette (ABC) transporters provide the first line of defense against xenobiotic chemicals and are upregulated in resistant pest organisms. While previous studies in bees have focused on the role of metabolic enzymes in insect detoxification, the presence and function of ABC transporters across the hive caste system remains largely unexplored. This study investigated the gene expression profiles of 12 ABC transporters known to be involved in chemical detoxification in arthropods across ten honey bee castes and life stages using quantitative real-time PCR. Protein homology to known MDR transporters in humans and Drosophila was inferred from BLAST and through phylogenetic analysis. Seven ABC genes that showed increased gene expression during worker bee development were identified as MDR-like transporters (AmeABCB1, AmeABCB6, AmeABCC1, AmeABCC4a-c, AmeABCG1); and their expression levels were further investigated in reproductive caste members (drone larvae, adult drones, queen ovaries, and queens). Significant variations were observed in defense gene expression among all castes suggesting reduced chemical defense capabilities in queens as evidenced by a dramatically reduced expression of five MDR-like transporter genes in queen bees relative to worker eggs: ABCB1 (4-fold), ABCC1 (2-fold), ABCC4a (2-fold), ABCC4b (3-fold), and ABCC4c (2-fold). Although our findings suggest that drones and queens are more vulnerable to direct xenobiotic exposure compared to workers, further research is required to better understand the different hive members’ responses to chemical threats.

Sample Acquisition

Honey bee samples were collected from hives maintained by the USDA-ARS Pollinator Health Research Unit in Davis, CA. Table S1 describes the sampling seasons, years, and locations. Qualitative monitoring of parasitic Varroa destructor mites was conducted using alcohol washes²¹, and levels were found to be low throughout the course of sample collection. Colonies were maintained according to standard best beekeeping practices in the United States²², and levels of Varroa destructor parasitic mites were maintained using treatments with Formic Pro^TM (NOD Apiary Products Ltd., Trenton, Ontario).

Ten different life stages were collected for this project:

·      Eggs (E)

·      Worker larvae instar 3 (WL3)

·      Worker larvae instar 5 (WL5)

·      Drone larvae instar 5 (DL5)

·      Worker pupae (P)

·      Newly emerged worker bees (NEB)

·      Nurses (N)

·      Foragers (F)

·      Drones (D)

·      Queens (Q)

Five samples of each life stage were collected for biological replication purposes. Additional queens were collected for dissection of queen ovaries.

Eggs (E) were collected by taking a frame containing E (from a hive or a queen monitoring cage²³ used in the laboratory) and tapping the E out onto a piece of dark-colored paper before pouring them into an RNase-free 1.5 mL tube. When E were sticking to the cells, they were removed using grafting tools. For each sample, 5-10 mg of E were collected.

All larvae and pupae were collected by grafting directly from hive frames and were stored in the same way. Adult bees (NEBs, nurses, foragers, and drones) were identified on live frames, collected, and euthanized by flash-freezing. During an annual queen replacement event in April 2022 mated queens were collected in the field and immediately flash-frozen in liquid nitrogen. All bee samples collected in the field were stored on dry ice until they could be brought back to the lab and stored at -80°C.

A subset of 5 queen samples from the previously collected mated queens had their ovaries isolated as an eleventh treatment group (QO) for the experiment. Dissections were performed using a Leica S9i stereo microscope (Leica Microsystems, Deerfield, IL). Frozen queens were dissected in a petri dish containing cool phosphate-buffered saline (PBS) (pH 7.4) solution on top of an ice block. 

Tissue Homogenization and RNA Extractions

Honey bee tissues were homogenized using the BeadBug™ 6 Six-Position Homogenizer (Benchmark Scientific, Sayreville, NJ, USA). Steel beads of 2.8 mm (Omni International, Kennesaw, Georgia, USA) were selected over ceramic beads to effectively break through the insects' exoskeletons and homogenize the samples. For each sample (except QO), the entire bee was used for RNA extraction. All samples were homogenized in DNA/RNA Protection Reagent (New England Biolabs, Ipswich, MA, USA) for three cycles of 30 seconds at a speed of 6 m/s with 45-second dwell periods. For samples with especially tough exoskeletons or wide bee samples (i.e., drones, queens) that impeded steel bead movement within the tube, up to three additional runs of the homogenization program were performed to ensure thorough soft tissue disruption. The sample was placed on ice between each run to prevent overheating that could potentially compromise RNA quality. RNA extraction was conducted using the Monarch® Total RNA Miniprep Kit (New England Biolabs, Ipswich, MA, USA). The quantities of reagents used in the initial steps varied based on sample weights. For instance, each WL3 sample weighing approximately 10 mg required minimal buffer and Proteinase K, whereas Q and D samples exceeding 200 mg necessitated maximal reagent usage. An optional on-column DNase 1 treatment was performed on every sample to ensure complete genomic DNA (gDNA) removal. RNA quality was validated using standard agarose gel electrophoresis and UV absorbance for determining A₂₆₀/A₂₃₀ ratios.

cDNA Conversion

Total RNA was converted into cDNA for qPCR analysis using the SuperScript IV First-Strand Synthesis System (Invitrogen, Thermo Fisher Scientific, Waltham, MA). The manufacturer's protocol was followed to prepare 1000 ng of cDNA per sample. If necessary, RNA samples were diluted with molecular biology-grade water before cDNA conversions to ensure accurate pipetting. For this reason, only 200 ng of cDNA was created for two of the E samples with particularly low RNA yield. Although the initial cDNA quantity differs between 1000 ng and 200 ng samples (while maintaining the same volume), the 1000 ng cDNA samples were diluted 1:4 to standardize the concentrations for all samples and to achieve the optimal concentration for qPCR.

Primer Design and Quality Control

The qPCR primers were designed using KEGG, CLC Genomic Workbench, Primer-BLAST (NCBI), and Primer3Plus²⁴. The gene sequences for all genes of interest, including their respective isoforms, were downloaded into CLC. Within each gene, isoforms were aligned, and gaps were recorded. Primers were subsequently designed within regions of the alignments where no gaps were present to capture all isoforms of each gene. Target amplicons were set between 80 and 200 base pairs (bp), with an optimal size of 120 bp. The potential for primer dimerization was assessed using Beacon Designer (Premier Biosoft International, Palo Alto, CA), and all primers were screened for potential non-specific binding using Nucleotide BLAST (NCBI). A previously published primer was used for the housekeeping gene RAD1A²⁵. For this project, three reference genes were selected: glyceraldehyde 3-phosphate dehydrogenase (GAPDH), ribosomal protein S5 (RS5), and Ras-related protein Rab-1A (RAD1A). The first two of these genes are commonly used as reference genes for honey bees, and the latter was recently identified as the most stable reference gene for qPCR experiments across honey bee life stages^25,26. High-purity salt-free (HPSF) primers (Eurofins Genomics, Louisville, Kentucky, U.S.) were reconstituted and diluted using TE buffer (pH = 7.4). The working primer concentrations were 10 µM.

Each primer pair specificity was assessed with nurse bee cDNA using standard PCR and agarose gel electrophoresis. Prior to qPCR analysis, the efficiency of the primers was also tested by creating a five-point cDNA serial dilution curve. The slope obtained from plotting concentrations against Ct values was used to calculate efficiency using the following equation²⁷:

[10^(-1/(slope))-1]*100

Efficiencies of the validated qPCR primers can be found in the supplemental table (Table S2).

Phylogeny

The protein sequence for each honey bee transporter gene was blasted using NCBI Protein BLAST against human and Drosophila sequences. The results were sorted by Max score, selecting only the top annotated sequence from each organism. Sequences for each organism were downloaded and organized in CLC Main Workbench 21.0.5 to construct a phylogenetic tree using the following parameters: Neighbor-Joining (algorithm), Jukes-Cantor (distance measure), and 1000 replicates (bootstrap). The root (outgroup) was set to be human HRas GTPase (NP_001123914).

Real-time – PCR (qPCR)

All qPCR plates were designed with three technical replicates for each unique reaction. Each plate contained one biological replicate of three different life stages, with five genes of interest and three reference genes. The 20 µL reaction protocol using SsoAdvanced Universal SYBR® Green Supermix (Bio-Rad, Hercules, CA, USA) was applied to all samples and pipetted onto MicroAmp™ EnduraPlate™ optical 96-well fast clear reaction plates (Thermo Fisher Scientific, Waltham, MA). After loading the reactions, the plates were sealed with MicroAmp™ optical adhesive film (Thermo Fisher Scientific, Waltham, MA) and centrifuged for one minute to thoroughly mix all the reagents.

All plates were run on a QuantStudio 3 Real-Time PCR system (Thermo Fisher Scientific, Waltham, MA) using the following protocol: a 30s hold at 98°C for polymerase activation and initial DNA denaturation (rate 4.3°C/s), followed by 40 cycles of 15s at 98°C for denaturation (rate 3°C/s), and 30s at 60°C for annealing/extension + plate read (rate 3°C/s). A melt-curve analysis was added after the 40 cycles: 15s at 95°C (rate 1.6°C/s), followed by 1min at 60°C (rate 1.6°C/s), and then 15s at 95°C (rate 0.05°C/s).

Using the biological replicates (n=5) ΔC_T values and technical plate replicates, an average value ΔC_T was obtained for each of the eleven treatments (E, WL3, WL5, DL5, P, NEB, N, F, D, Q, QO) and twelve genes of interest (Table 1). A geometric mean of all housekeeping genes was used to calculate the ΔC_T values. The qPCR data was analyzed using the Livak method of calculating fold change with ΔΔC_T ²⁸, comparing differences between treatments. Each tested life stage was considered a "treatment" except for E, which was designated as the "control" treatment. E was selected as the control group since all caste members pass through this stage. Gene expression was normalized to its corresponding expression in E (= value of 1). The fold change of gene expression was calculated using the equation: 

2^e- [ΔΔCT = ΔCT ("treatment") – ΔCT ("control")]

All qPCR results were assessed for quality criteria, including melting temperature, C_T values greater than 35, and technical replicate C_Ts within 0.5 units of their mean. Quality parameters were applied using RStudio v.4.2.2 to ensure unbiased analysis. The filtered data were then used for statistical and interpretation.

Statistical analysis

All the statistical analysis was performed using RStudio v.4.2.2. The filtered data were first checked for normality using a histogram, QQ-plot, and a Shapiro-Wilk normality test (W=0.97, p-value < 2.2 e^-16). The normality test and the models created to analyze the data were applied to the housekeeping normalized C_T values (ΔC_T). A one-way  ANOVA test (“analysis of variance”) was conducted to determine if significant differences were present across the twelve treatments (E, WL3, WL5, DL5, P, NEB, N, F, D, Q, QO) for each gene of interest. Following the ANOVA, a Tukey test or Tukey’s Honest Significant Difference test was applied to the data to identify significant differences between combinations of treatments. The confidence interval for the Tukey model was set to 0.99. Comparisons were made across the same gene and across the same life stage, with values below 0.01 considered statistically significant.