DOI: 10.5061/dryad.d2547d8c2
title: Supplemental information and raw data: Developmental and caste-specific expression patterns of ATP-Binding Cassette (ABC) transporters in honey bees (Apis mellifera)

description: 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 &gt; drones &gt; queens) and
 suggest caste-specific trade-offs between reproduction and chemical
 defense in honey bee superorganisms.  
description: 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 <em>Varroa
 destructor</em> 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 <em>Varroa
 destructor</em> 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
 A<sub>260</sub>/A<sub>230</sub> and
 A<sub>260</sub>/A<sub>280</sub> 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 <em>Drosophila melanogaster</em> . Full-length
 protein sequences of honey bee ABC transporters were used to identify
 corresponding orthologs in human and <em>Drosophila</em>
 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) ΔC<sub>T</sub>
 values and technical plate replicates, an average
 ΔC<sub>T</sub> 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 ΔC<sub>T</sub>
 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.
 <strong>Protein
 name</strong> <strong>BeeBase/HGD
 ID</strong> <strong>Gene ID</strong>
 <strong>Gene
 symbol</strong> <strong>Gene
 description</strong> <strong>Length
 (aa)</strong> 
 
 <em>Ame</em>-ABCB1 GB55378 551167 Mdr49 Multidrug resistance protein homolog
 49 1343 
 
 <em>Ame</em>-ABCB6 GB54214 726867 Hmt−1 ABC transporter ATP-binding protein/permease
 Hmt−1 840
 
 <em>Ame</em>-ABCC1 GB53134 412217 MRP Multidrug-Resistance like Protein 1
 1536 
 
 <em>Ame</em>-ABCC10 GB44193 725993 LOC725993 Multidrug resistance-associated protein
 7 1625 
 
 <em>Ame</em>-ABCC4a GB50195 413959 LOC413959 Multidrug resistance-associated protein
 4 1356 
 *Ame-*ABCC4b GB45277 551061 LOC551061 Multidrug resistance-associated protein
 4 1392 
 
 <em>Ame</em>-ABCC4c GB45278 725051 LOC725051 Multidrug resistance-associated protein
 4 1418 
 
 <em>Ame</em>-ABCC4d GB17987; GB43189 410269 LOC410269 Probable multidrug resistance-associated protein
 lethal (2) 03659 1333 <em>Ame</em>-ABCG1
 GB44770 412748 LOC412748 ATP-binding cassette sub-family G member
 1 631 
 
 <em>Ame</em>-ABCG2a GB49098 410967 E23 Early gene at 23 657 *Ame-*ABCG2b GB43076 726508 LOC726508 Protein scarlet 623 <em>Ame</em>-ABCG5
 GB41827 726513 LOC726513 ATP-binding cassette sub-family G member
 5 625 
 All qPCR results were assessed for
 quality criteria, including melting temperature, C<sub>T</sub>
 values greater than 35, and technical replicate C<sub>T</sub>s
 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
 &lt; 2.2 e^−16^). Normality testing and models created to analyze the
 data were applied to housekeeping gene-normalized
 C<sub>T</sub> values (ΔC<sub>T</sub>). 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&lt;0.05 or p&lt;0.05, p&lt;0.01, p&lt;0.001, and
 non-significance (ns).
description: # Supplemental information and raw data: Developmental and caste-specific
 expression patterns of ATP-Binding Cassette (ABC) transporters in honey
 bees (Apis mellifera)
 [https://doi.org/10.5061/dryad.d2547d8c2](https://doi.org/10.5061/dryad.d2547d8c2) ## Description of the data and file structure Supplemental Information and qPCR raw data from: Developmental and Caste-specific Expression Patterns of ATP-Binding Cassette (ABC) Transporters in Honey bees (*Apis mellifera*) ### Files and variables #### Files: #### LS_qPCR_Organizeddata_All.csv \**Description: This file contains the raw c for all ABC transporters, and housekeeping genes. ##### Variables and descriptions * null = no data was reported by the qPCR instrument for this cell.  - sample_name (categorical) = name of the sample specimen and biological replicate number, for example D1 is drone biological replicate 1, while D2 would be drone biological replicate 2.  * D = Drone bee * F = Forager bee  * N = Nurse bee  * NEB = Newly emerged bee  * Q = Queen bee * QO = Queen bee ovaries * QL = Queen bee Larvae * QP = Queen bee Pupa * E = Honey bee worker eggs * L3 = Worker bee Instar larva 3 * L5 = Worker bee Instar larva 5 * P = Worker bee pupa * DL5 = Drone bee Instar larva 5 - target_genes (categorical) = name of the gene of interest  * ABC = ATP biding cassette transporter * GAPDH =Glyceraldehyde 3-phosphate dehydrogenase * RAD1 =Cell cycle checkpoint protein RAD1 * RS5 =Ribosomal Protein S5 * Scarlet =pigment transporter gene - ct_values (numeric) = numeric value obtained from qPCR analysis  - Tm (numeric) = Melting temperature values obtained from melt curve analysis in qPCR instrument. Each sample should only have one melting temperature value (Tm1). If the sample has more than one temperature value (Tm2) then that sample is excluded from further statistical analysis due to possible contamination and low reliability of that ct value.  #### LS_FilteredData.csv \**Description: This file contains the filtered qPCR raw data. The data was filtered based on three exclusion parameters: 1. Samples showing a Tm2 value. 2. Samples with a ct value higher or equal to 35. 3. Samples with a ct values of its technical replicates which were more than 0.5 units apart from the mean of those replicates. ##### Variables and descriptions * sample_name (categorical) = name of the sample specimen and biological replicate number, for example D1 is drone biological replicate 1, while D2 would be drone biological replicate 2.  * D = Drone bee * F = Forager bee  * N = Nurse bee  * NEB = Newly emerged bee  * Q = Queen bee * QO = Queen bee Ovaries * QL = Queen bee Larvae * QP = Queen bee Pupa * E = Honey bee worker eggs * L3 = Worker bee Instar larva 3 * L5 = Worker bee Instar larva 5 * P = Worker bee pupa * DL5 = Drone bee Instar larva 5 * target_genes (categorical) = name of the gene of interest  * ABC = ATP biding cassette transporter * GAPDH =Glyceraldehyde 3-phosphate dehydrogenase * RAD1 =Cell cycle checkpoint protein RAD1 * RS5 =Ribosomal Protein S5 * Scarlet =pigment transporter gene * ct_values (numeric) = numeric value obtained from qPCR analysis  * Tm (numeric) = Melting temperature values obtained from melt curve analysis in qPCR instrument. Each sample should only have one melting temperature value (Tm1). * n (numeric) = number of technical replicates.  * Mean (numeric) = Average value of the technical replicates.  * Stdev (numeric) = Standard deviation value of the technical replicates.  * RSD (numeric) = Relative standard deviation value of the technical replicates.  #### Supplemental_Bee_ABCTransporters.docx \**Description: This file contains the supplemental figures and tables for the full publication ([https://doi.org/10.1016/j.etap.2025.104789)](https://doi.org/10.1016/j.etap.2025.104789 "Persistent link using digital object identifier") ## Code/software Microsoft Excel or OpenOffice or LibreOffice ## Access information Other publicly accessible locations of the data: * NA Data was derived from the following sources: * NA ## Version changes **04-sep-2025:** Added additional tables to the methods section (table 1, table 2). Added two new lifestages for queen bee and its respective data set (QP, QL). Added the supplemental data file (Supplemental_Bee_ABCTransporters.docx).