Supplemental information and raw data: Developmental and caste-specific expression patterns of ATP-Binding Cassette (ABC) transporters in honey bees (Apis mellifera)

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. This study investigated the gene expression profiles of 12 ABC transporters involved in chemical detoxification across three honey bee castes and 13 life stages using quantitative real-time PCR. Six ABC genes showed increased expression during worker bee development and were identified as MDR-like transporters (Ame-ABCB1, Ame-ABCB6, Ame-ABCC4a-c, Ame-ABCG1). Four transporters showed pupal-specific expression during metamorphosis. Queens exhibited dramatically reduced MDR transporter expression compared to workers between 1.7-fold lower (ABCB6) and 17.5-fold lower (ABCB1). Drones showed intermediate expression levels. Queen ovaries demonstrated tissue-specific upregulation of select transporters. These findings reveal a vulnerability hierarchy (foragers > drones > queens) and suggest caste-specific trade-offs between reproduction and chemical defense in honey bee superorganisms.

Methods
2.1. Sample acquisition

Honey bee samples were collected from hives maintained by the USDA-ARS Pollinator Health Research Unit in Davis, CA. Table S2 describes the sampling seasons, years, and locations. Quantitative monitoring of parasitic Varroa destructor mites was conducted using alcohol washes (Jack and Ellis, 2021), and levels were found to be low throughout sample collection. Colonies were maintained according to standard best beekeeping practices in the United States (Best Management Practices for Beekeepers and Growers – Bee Health, 2024), and levels of Varroa destructor mites were controlled using treatments with Formic Pro™ (NOD Apiary Products Ltd., Trenton, Ontario).

Five samples of each life stage were collected for biological replication. Additional queens were collected for dissection of five ovaries. Eggs (E) were collected by taking a frame containing E (from a hive or queen monitoring cage (Fine et al., 2021) used in the laboratory) and tapping them out onto 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 similarly.

Adult bees (NEBs, nurses, foragers, and drones) were identified on live frames, collected, and euthanized by flash-freezing in liquid nitrogen. During an annual queen replacement event in April 2022, mated queens were collected in the field and immediately put on dry ice. All bee samples were stored at −80°C until processing. A subset of five queen samples from these mated queens had their ovaries isolated as an additional treatment group (QO). All dissections were performed using a Leica S9i stereo microscope (Leica Microsystems, Deerfield, IL). Frozen queens were dissected in a petri dish containing chilled phosphate-buffered saline (PBS) (pH 7.4) over ice.

2.2. Tissue homogenization and RNA extraction

Honey bee tissues were homogenized using BeadBug™ 6 Six-Position Homogenizer (Benchmark Scientific, Sayreville, NJ). Steel beads of 2.8 mm (Omni International, Kennesaw, Georgia) were selected over ceramic beads to effectively break through insect exoskeletons and homogenize the samples. For each sample (except QO), the entire bee was used for total RNA extraction. All samples were homogenized in DNA/RNA Protection Reagent (New England Biolabs, Ipswich, MA) for three cycles of 30 s at 6 m/s with 45-second dwell periods. For samples with especially tough exoskeletons or large specimens (i.e., drones and queens) that impeded steel bead movement within the tube, up to three additional homogenization cycles were performed to ensure thorough soft tissue disruption. Samples were placed on ice between each cycle 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). The quantities of reagents used in the initial steps varied based on sample weights. For example, each WL3 sample (weighing approximately 10 mg) required minimal buffer and Proteinase K, whereas Q and D samples (exceeding 200 mg) necessitated maximum reagent usage. The optional on-column DNase I treatment was performed on every sample to ensure complete genomic DNA (gDNA) removal. RNA quality was validated using standard 2 % w/v agarose gel electrophoresis and UV absorbance to determine the A260/A230 and A260/A280 ratios (Spooner et al., n.d).

2.3. cDNA conversion

Total RNA was converted to 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, except for E samples, which had particularly low RNA yields (200 ng). For qPCR analysis, all final cDNA concentrations were standardized to 20 ng/µl.

2.4. Primer design and quality control

All qPCR primers were designed using KEGG, CLC Genomic Workbench v20.0.4, Primer-BLAST (NCBI), and Primer3Plus (Spooner et al., n.d). 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. 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). For this project, three reference genes were selected: glyceraldehyde 3-phosphate dehydrogenase (GAPDH), ribosomal protein S5 (RS5), and Ras-related protein Rab-1A (RAB1A). The first two 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 (Kim et al., 2021, Kim et al., 2022). The previously published primer pair was used for the housekeeping gene RAB1A (Kim et al., 2021). High-purity salt-free (HPSF) primers (Eurofins Genomics, Louisville, Kentucky) were reconstituted and diluted using TE buffer (pH = 7.4). Working primer concentrations were 10 µM.

The specificity of each primer pair was assessed using nurse bee cDNA via standard PCR and agarose gel electrophoresis (Spooner et al., n.d). Prior to qPCR analysis, primer efficiency was tested and optimized 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 (Rogers-Broadway and Karteris, 2015):

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

Primer pairs with reaction efficiencies between 90 % and 110 % were considered acceptable. Calculated efficiencies of all validated qPCR primers can be found in Table S3.

2.5. Phylogeny

Phylogenetic analysis was performed to determine the evolutionary relationships of honey bee ABC transporters with their MDR-like orthologs in human and Drosophila melanogaster . Full-length protein sequences of honey bee ABC transporters were used to identify corresponding orthologs in human and Drosophila through reciprocal NCBI BLAST. The BLAST 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 v21.0.5 to perform multiple sequence alignments and construct a phylogenetic tree using the following parameters: Neighbor-Joining (algorithm), Jukes-Cantor (distance measure), and 1000 replicates (bootstrap). The tree was rooted with the human HRas GTPase (NP_001123914) as an outgroup.

2.6. Real-time PCR (qPCR)

All qPCR plates were designed with three technical replicates for each unique reaction. Each plate contained one biological replicate from three different life stages, with five genes of interest and three reference genes. A 20 µL reaction protocol using SsoAdvanced Universal SYBR® Green Supermix (Bio-Rad, Hercules, CA) 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, plates were sealed with MicroAmp™ optical adhesive film (Thermo Fisher Scientific, Waltham, MA) and centrifuged for one minute to thoroughly mix all reagents.

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

Using biological replicates (n = 5) ΔCT values and technical plate replicates, an average ΔCT was obtained for each of the eleven treatments (Table 1) and twelve genes of interest (Table 2). A geometric mean of all three housekeeping genes was used to calculate ΔCT values. qPCR data were analyzed using the Livak method for fold change calculation using ΔΔC~T ~(Spooner et al., n.d; Livak and Schmittgen, 2001), comparing differences between treatments. Each tested caste member, developmental stage, and tissue type were considered a "treatment" except for E, which was designated as the "control" since all caste members pass through this stage. Gene expression was normalized to its corresponding expression in E (= value of 1). Fold change was calculated using the following equation:

Table 1. Honey bee samples collected for qPCR analysis of ABC transporter gene expression by caste, developmental stage, and tissue type. The sample abbreviations are used consistently throughout the text. Eggs represent those destined to develop into workers. For details on adult bee sample collection, please refer to the Supplemental Table S1.

Caste
Life Stage/Tissue
Abbreviation

Worker bee
Eggs
E

Larvae (Instar 3)
WL3

Larvae (Instar 5)
WL5

Pupae
P

Newly Emerged Bee
NEB

Nurse
N

Forager
F

Drone bee
Larvae (Instar 5)
DL5

Adult
D

Queen bee
Ovaries
QO

Larvae
QL

Pupae
QP

Adult
Q

Table 2. Metadata summary of the 12 honey bee ABC transporters investigated in this study. Gene ID and naming conventions in NCBI and BeeBase/HGD for honey bee ABC transporters are shown. The four different genes for honey bee ABCC4 (a-d) are all paralogous to human ABCC4. Two ABCG (a-b) genes were paralogous to human ABCG2.

Protein name
BeeBase/HGD ID
Gene ID
Gene symbol
Gene description
Length (aa)

Ame-ABCB1
GB55378
551167
Mdr49
Multidrug resistance protein homolog 49
1343

Ame-ABCB6
GB54214
726867
Hmt−1
ABC transporter ATP-binding protein/permease Hmt−1
840

Ame-ABCC1
GB53134
412217
MRP
Multidrug-Resistance like Protein 1
1536

Ame-ABCC10
GB44193
725993
LOC725993
Multidrug resistance-associated protein 7
1625

Ame-ABCC4a
GB50195
413959
LOC413959
Multidrug resistance-associated protein 4
1356

*Ame-*ABCC4b
GB45277
551061
LOC551061
Multidrug resistance-associated protein 4
1392

Ame-ABCC4c
GB45278
725051
LOC725051
Multidrug resistance-associated protein 4
1418

Ame-ABCC4d
GB17987; GB43189
410269
LOC410269
Probable multidrug resistance-associated protein lethal (2) 03659
1333

Ame-ABCG1
GB44770
412748
LOC412748
ATP-binding cassette sub-family G member 1
631

Ame-ABCG2a
GB49098
410967
E23
Early gene at 23
657

*Ame-*ABCG2b
GB43076
726508
LOC726508
Protein scarlet
623

Ame-ABCG5
GB41827
726513
LOC726513
ATP-binding cassette sub-family G member 5
625

All qPCR results were assessed for quality criteria, including melting temperature, CT values greater than 35, and technical replicate CTs within 0.5 units of their mean. Quality parameters were applied using RStudio v.4.2.2 (R-Packages: dplyr, tidyverse, openxlsx) to ensure unbiased analysis. Filtered data were then used for statistical interpretation.

2.7. Statistical analysis

All statistical analysis was performed using RStudio v.4.2.2 (R-packages: dplyr, tibble, forcats, emmeans, multcomp, multcompView, ggplot2, openxlsx). Filtered data were first checked for normality using histogram, QQ-plot, and Shapiro-Wilk normality test (W=0.97, p-value < 2.2 e^−16^). Normality testing and models created to analyze the data were applied to housekeeping gene-normalized CT values (ΔCT). A one-way ANOVA (“analysis of variance”) was conducted to determine whether significant differences were present across the thirteen treatments (E, WL3, WL5, DL5, P, NEB, N, F, D, Q, QO, Ql, QP) for each gene of interest. Following ANOVA, a post-hoc pairwise comparison between sample groups were performed using least-squares means with significance levels of p<0.05 or p<0.05, p<0.01, p<0.001, and non-significance (ns).