Acute Sublethal Effects of the Neonicotinoid Imidacloprid
on the Honeybee Brain Transcriptome

A THESIS SUBMITTED TO THE SCHOOL OF GRADUATE STUDIES OF
BLOOMSBURG UNIVERSITY OF PENNSYLVANIA

IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
MASTERS OF BIOLOGY

DEPARTMENT OF BIOLOGICAL AND ALLIED HEALTH SCIENCES

BY
HEATHER J. LLEWELLYN

BLOOMSBURG, PENNSYLVANIA
SUMMER 2020

Abstract
Global declines in honeybees have been linked to widespread use of pesticides.
Sublethal doses of the neonicotinoid pesticide imidacloprid have been shown to cause
physiological and behavioral changes that negatively impact hive health. This study
examined the effects of acute, sublethal doses of imidacloprid on the honeybee gene
expression in the brain.
Experiment 1 identified an imidacloprid dosage that yielded a cellular stress
response in honeybees. Honeybee foragers were harnessed, fed to satiation, and randomly
assigned to control (1.5 M sucrose) or treatment groups receiving imidacloprid at doses
of 1/5th, 1/10th, 1/20th, 1/50th, 1/100th, and 1/500th of the LD50 (18.0 ng/bee). After four
hours, impaired motor responses and elevated levels of Heat Shock Protein 70 and
Superoxide Dismutase, markers of cellular and oxidative stress, respectively, were
observed at sublethal imidacloprid doses. A conservative dose of 0.9 ng/bee (1/20th of the
LD50) was selected to treat bees for RNA transcriptome analysis in Experiment 2.
In Experiment 2, bees were randomly assigned to four acute-exposure groups:
Control-0h, Control-4h, Treatment-0h, and Treatment-4h. Control bees were fed 1.5 M
sucrose while treatment groups received 0.9 ng/bee of imidacloprid. RNA isolated from
bee brain tissue was sent to the University of Illinois at Urbana-Champaign for RNA
sequencing. Of the 10,597 genes recovered from the reference genome, 4,205 genes were
differentially expressed. Collectively, the Differential Expression Analysis,
Multidimensional Scaling Plot, and heat map agree that the Control-4h group had the
greatest changes in gene expression, counter to our prediction that the greatest alteration
would be in imidacloprid-treated bees. Corroborating evidence supported the post-hoc
hypothesis that the Control-4h and the Treatment-4h samples were switched and
mislabeled prior to shipping. Comparisons between Control-4h and Treatment 4-h groups
remain valid. If samples were switched, only the direction of differential gene expression
(i.e., up-regulation versus down-regulation) would be affected by imidacloprid. Gene set
enrichment analysis indicated that the key pathways affected were: oxidative
phosphorylation, longevity regulating pathway, apoptosis, peroxisome, FOXO signaling,
drug metabolism- cytochrome P450, metabolism of xenobiotics by cytochrome P450,
circadian rhythm, and glutathione metabolism. These gene networks relate to key
2

biological functions of honeybees that have the potential to affect colony viability. Future
research will focus on hypothesis-driven gene expression studies that relate specific
molecular changes to biological functions and organism-level performance, an integrative
approach that is essential to understanding the declines of these essential pollinators.

3

Acknowledgments
After four years, the journey has finally come to an end. This research would not
have been possible without Aucker’s Apiary in Millville, Pennsylvania and Dr. John M.
Hranitz for allowing the collection of honeybees from their hives. Joshua I. Petersheim
spent countless hours collecting, harnessing and feeding, treating, and capturing escaped
honeybees. Erin Smith worked towards the development and validation of biochemical
assays for catalase and superoxide dismutase to compare oxidative stress in control and
imidacloprid-treated honeybees. Both Joshua and Erin held integral roles towards the
completion of this research.
The Bloomsburg University Professional Experience Grant (2019 and 2020),
Bloomsburg University Graduate School, and the Charlotte P. Magnum Student Support
Program allowed for travel and board to the Society for Integrative and Comparative
Biology National Conference in Tampa, Florida and Austin, Texas. If it was not for this
funding, priceless opportunities to present and listen to research would have been missed.
The Research and Scholarship Grant from Bloomsburg University allowed for the
purchase of necessary supplies to complete this thesis research. The Department of
Biological and Allied Health Sciences and College of Science and Technology at
Bloomsburg University of Pennsylvania which allowed me to be a part of this academic
community, I am forever grateful.
My Thesis Committee: Drs. Surmacz, Hare-Harris, Hranitz, and Schwindinger for
their time, patience, and willingness towards the completion of this research. My
schedule was not easy to work with but yet, at any day of the week, at all hours, we met
to discuss and work on all the moving parts of this project. The Keck Center for
Comparative and Functional Genomics at the University of Illinois Urbana-Champaign
for graciously rerunning our samples for RNA sequencing. Dr. Hare-Harris for
painstakingly writing a python and R script, data cleaning, and reviewing DAVID/KEGG
analysis. If it was not for you, I would still be trying to open a file and figure out where I
sent my python output. Dr. Hranitz who taught me more than I ever thought I would need
to know about honeybees. Dr. Klinger for calling me that fateful day in the spring of
2016 to let me know there was an advanced DNA class at 2 p.m. and that I should attend
because I was accepted into the Graduate program; thank you for giving me the chance.
5

Dr. Surmacz, I honestly do not know what to say about this marvelous educator. She has
been present for all of the highs and lows, laughs and tears, breakdowns and successes of
this research. These four years have been a blessing for me. You have shown me how to
persevere through any type of situation with a smile and kind words or just a simple
“buck up buttercup, that’s research.” You have made a huge impact.
A huge thank you to my family and friends, you have allowed me to be a
complete nut-case for four years. Mom, you are and will always be my biggest fan.
Thank you for always believing I can do anything and everything even if I failed. My
future husband, Zachary, thank you for always being there, wishing me luck, eating my
feelings with me, and telling me to get it together even when I could not. Lastly, to all
those who said it could not be done. It can and it was.

6

Table of Contents
Abstract…………………………………………………………………………….. 2
Acknowledgments…………………………………………………………………..5
Table of Contents…………………………………………………………………... 7
List of Figures ……………………………………………………………………... 10
List of Tables………………………………………………………………………..12
List of Appendices…………………………………………………………………. 13
Introduction………………………………………………………………………… 14-38
Biology of the Honeybee…………………………………………………... 14
The Honeybee Genome……………………………………………………..18-24
Genome Organization ……………………………………………... 19
Functional Categories of Genes……………………………………. 20-24
Immune System Pathways…………………………………. 20
Anti-Xenobiotic Defense Mechanisms…………………….21
Antioxidants………………………………………………... 22
Heat Shock Protein………………………………………… 22
Circadian Rhythm and Behavior……………………………23
Colony Collapse Disorder………………………………………………….. 24
Pesticides as Environmental Stressors……………………………………... 25
Pesticide Dynamics in Honeybees…………………………………………. 29
Effects of Neonicotinoid Pesticides on the Honeybee Transcriptome……...31
Sublethal Stress in Honeybees……………………………………………... 32
Cellular Responses to Stressors by Honeybees……………………………. 35
Oxidative Stress……………………………………………………………. 37
Objectives and Research Questions………………………………………………... 38
7

Methods……………………………………………………………………………. 40-58
Experiment 1……………………………………………………………….. 40-47
Experimental Design……………………………………………….. 40
Honeybee Collection………………………………………………. 41
Motor Testing……………………………………………………….43
Homogenization of Bee Head Capsules…………………………… 43
Superoxide Dismutase Assay………………………………………. 43
Protein Assay………………………………………………………. 45
HSP70 ELISA ……………………………………………………... 45
Statistics……………………………………………………………. 47
Experiment 2……………………………………………………………….. 47-58
Experimental Design………………………………………………. 47
Honeybee Collection and Treatments ……………………………... 48
Head Capsule Dissection……………………………………………48
RNA Isolation ………………………………………………………49
RNA Sequencing……………………………………………………53
Post-Sequencing Analysis and Statistics……………………………54
Gene Enrichment and Pathway Analysis…………………………... 56
Results ………………………………………………………………………………58-68
Experiment 1……………………………………………………………….. 58-62
Motor Function Response of Honeybees…………………………... 58
Superoxide Dismutase Activity in Honeybees…………………….. 60
HSP70 Concentrations in Honeybees ………………………………61
Overall Results of Experiment 1……………………………………62

8

Experiment 2……………………………………………………………….. 62-68
Gene Expression Analysis: RNA-Sequencing……………………... 62
Discussion………………………………………………………………………….. 68-83
Motor Function Responses………………………………………………… 68
Superoxide Dismutase Activity……………………………………………. 70
HSP70 ………………………………………………………………………71
Determination of Imidacloprid Dose to Use for Experiment 2……………..72
Mortality…………………………………………………………………… 72
RNA-Seq Results…………………………………………………………... 72
Post-Hoc Hypothesis and Supporting Evidence…………………………… 73
Pathway and Gene Set Enrichment Functional Analysis…………………... 75
Limitations…………………………………………………………………………. 83
Future Directions…………………………………………………………………... 85
Conclusion…………………………………………………………………………. 86
Literature Cited…………………………………………………………………….. 89
Appendices…………………………………………………………………………. 100

9

List of Figures
Figure 1: Differentiation of body type between classes of bees within a hive…….14
Figure 2: An adult honeybee……………………………………………………….15
Figure 3: Life cycle of honeybees………………………………………………….17
Figure 4: Ideogram and karyotype for Apis mellifera………………………………….19
Figure 5: U.S. managed honeybee colony loss estimates………………………….24
Figure 6: Neurotransmission with acetylcholine and neonicotinoids………………27
Figure 7: Chemical structure and molecular formula of imidacloprid……………..29
Figure 8: Proposed mechanisms underlying the honeybee response to nicotine
exposure……………………………………………………………………………..31
Figure 9: Dose-response relationship for acute toxicity……………………………33
Figure 10: Sublethal Stress (SLS) Model…………………………………………..34
Figure 11: Hormesis dose-response curve………………………………………….36
Figure 12: HSP70 concentrations among pretreatment groups of honeybees and those
treated with ethanol…………………………………………………………………37
Figure 13: Experiment 1 design……………………………………………………41
Figure 14: Honeybee feeder, collection, and harnessing…………………………. 42
Figure 15: Honeybee treatment…………………………………………………… 43
Figure 16: Superoxide dismutase reaction and assay principle…………………... 44
Figure 17: Superoxide dismutase assay microplate………………………………. 44
Figure 18: Bio-Rad DC protein assay…………………………………………….. 45
Figure 19: Monoclonal antibody ELISA test for HSP70…………………………. 46
Figure 20: Experiment 2 design……………………………………………………48
Figure 21: Procedure for isolating the honeybee brain…………………………… 49
Figure 22: Schematic of RNA isolation procedure………………………………. .50
Figure 23: Formaldehyde gel for RNA…………………………………………… 52
10

Figure 24: TruSeq stranded mRNA library preparation workflow………………….53
Figure 25: Pairwise comparisons of control and imidacloprid-treated bees………. 56
Figure 26: Average motor scores of honeybees exposed to sublethal doses of
imidacloprid……………………………………………………………………….....59
Figure 27: Imidacloprid motor response curve……………………………………. 60
Figure 28: Superoxide dismutase activity in honeybees treated with
imidacloprid………………………………………………………………………....61
Figure 29: The effects of sublethal doses of imidacloprid on HSP70 levels in
honeybees……………………………………………………………………………62
Figure 30: One-way ANOVA heat map of Control-0h, Control-4h, Treatment-0h, and
Treatment-4h…………………………………………………………………………65
Figure 31: Multi-dimensional scaling plot performed on 5,000 of the most variable
genes…………………………………………………………………………………66
Figure 32: Peroxisome Pathway of Treatment-4h-Control-0h………………………76
Figure 33: Circadian Rhythm Pathway-Fly of Treatment-4h-Control-0h…………...77
Figure 34: FOXO Signaling Pathway of Treatment-4h-Control-0h…………………78

11

List of Tables
Table 1: Average absorption ratio and average protein present……………………52
Table 2: Number of genes removed for each comparison………………………… 56
Table 3: Differential Gene Expression (DE) for all four pairwise comparisons….. 64
Table 4: Key pathways and genes at Control-4h and Control-0h…………………. 67
Table 5: Removal of genes and outliers using Python…………………………….. 74
Table 6: Results of further data cleaning using Python…………………………….75
Table 7: Enriched genes at Treatment-4h and Control-0h………………………… 81

12

List of Appendices
Appendix A: Animal Research (IACUC) Approval……………………………………… 100
Appendix B: Literature data on neonicotinoid residues in bee-collected pollen, honey, and
bees………………………………………………………………………………………… 101
Appendix C: Lethal and sub-lethal effects by imidacloprid to individual (organism level) honey
bees as determined in different studies by oral exposure under laboratory conditions……. 102
Appendix D: Concentrations of imidacloprid causing lethal and sublethal effects on (micro)colony level in honey bees as determined in different studies by oral exposure under laboratory
conditions……………………………………………………………………………………103
Appendix E: Macho® Stock solution preparation…………………………………………104
Appendix F: Treatment group preparation…………………………………………………105
Appendix G: Bovine stock albumin standard preparation…………………………………106
Appendix H: Stock bovine HSP70 standard preparation………………………………… .107
Appendix I: Information needed to submit RNA for library construction and sequencing to
University of Illinois Keck Center………………………………………………………… 108
Appendix J: RNA-Seq Report prepared by the University of Illinois……………………. 111
Appendix K: Sequence quality from Multi-QC report……………………………………. 117
Appendix L: Outlier genes removed from each comparison group………………………..118
Appendix M: Python script……………………………………………………………….. 121
Appendix N: Average Motor Scores of honeybees collected and treated with
imidacloprid………………………………………………………………………………....122
Appendix O: List of 100 most significant changes in Treatment-4h versus Control-0h…..123
Appendix P: Control Time 4 Hour-Treatment Time 0 Hour KEGG pathways……………125
Appendix Q: Treatment Time 4 Hour-Control Time 0 Hour KEGG pathways……………135
Appendix R: DAVID Analysis of significant (p<0.1) enriched genes for Treatment-4h versus
Control-0h…………………………………………………………………………………...146
Appendix S: Enriched genes at Control-4h-Treatment-0h…………………………………176

13

I.

Introduction
Apis mellifera L., commonly known as the honeybee, is one of the world’s most

important pollinators both economically and ecologically. With honeybees pollinating 71
out of the 100 most common crops and accounting for 90% of the worlds food supply, it
is important to determine why honeybees and other pollinators are dying off at alarming
rates (Pilatic 2012, du Rand et al. 2015). Pesticides, malnutrition, habitat loss, parasites
and pathogens have all been identified as interacting factors that contribute to the drastic
loss of honeybees (du Rand et al. 2015).
The Biology of the Honeybee:
A bee colony contains two sexes, male and female (Figure 1). There are two
classes of female bees: queen bee and worker bees. Eggs laid by the queen develop into
mature honeybees. Unfertilized eggs develop into male drones (haploid), while fertilized
eggs develop into either female workers or queen cells (diploid). This type of sex
determination, haplodiploidism, is a characteristic of insects in the order Hymenoptera.
The metamorphosis from honeybee larva to pupa takes place within sealed cells and takes
16 days after deposition of the egg for the adult queen, 21 days for worker bees, and 24
days for drones (Agriculture and Consumer Protection 1990).

Figure 1: Differentiation of body type between classes of bees within a hive (Ellis and Mortenson
2017).

14

An adult queen is the reproductive female in the honeybee colony. Compared to a
worker bee, the queen has a longer and plumper abdomen with a head and thorax of
similar size (Ellis and Mortenson 2017). Worker honeybees possess a smaller body and
are specialized for pollen and nectar collection by having longer tongues and larger crops.
Worker bees have reduced ovaries, a stinger, and are not capable of mating. To assist in
carrying the large amounts of pollen back to the hive, the worker bee has a unique feature
called a corbicula or more commonly known as pollen basket (Ellis and Mortenson
2017). Drones, the only male caste within the hive, are easily differentiated from the
females. An adult drone possesses a larger thorax, no stinger and has “fly-like” eyes that
touch at the center at the top of their head. The drone’s abdomen is considered bulletshaped since it is thick with a blunt end (Ellis and Mortenson 2017). A general depiction
of an adult honeybee is seen in Figure 2 below.

Figure 2: An adult honeybee (Photo Credit: BuzzAboutBees 2010-2020). Covered in branched hairs, an
adult honeybee has three body regions: head, thorax and abdomen. The head contains compound eyes and
antennae. Two pairs of wings and three pairs of legs are attached to the thorax. In females, a barbed stinger
contains a poison-sac (Ellis and Mortenson 2017).

Each colony contains one fertile queen who lays eggs and 20,000-80,000 sterile
worker bees that maintain the colony. In comparison to the females, there are
approximately 300-800 fertile males, making the colony predominantly female
(Agriculture and Consumer Protection 1990). Among the mature honeybees there are
15

several thousand immature bees or brood; approximately 5,000 eggs and anywhere from
25-30,000 immature bees that are in various stages of development. Approximately
10,000 of the immature bees are newly hatched larvae (Agriculture and Consumer
Protection 1990).
A honeybee’s role within the hive is dependent on the food source fed to them. In
the first three days, all developing eggs are fed with “bee milk” or “royal jelly” produced
by nurse bees, younger worker bees that are not ready to leave the hive (Agriculture and
Consumer Protection 1990). After the first three days, worker and drone larvae are fed
mixed food that is composed of honey and pollen. Unlike the worker and drone larvae,
the queen is fed royal jelly for her entire larval life of approximately five days.
A honeybee egg measures between 1 to 1.5 mm in length and is often described as
appearing like a tiny grain of rice. These eggs can be found in individual hexagonal
shaped wax cells located in the brood area of the comb. The honeybee eggs hatch after
three days and larvae appear (Ellis and Mortenson 2017; Figure 3). Larvae are white in
color and appear as a curled “C” shape at the bottom of their wax cell. The amount of
time a honeybee remains in the larval stage is dependent on caste (worker: 6 days, drone:
6.5 days, and a queen: 5.5 days). As a pupa, the honeybee body becomes extended into an
upright position within the cell which is cover by the adult workers with a wax cap (Ellis
and Mortenson 2017). Pupa remain under the wax cap until it is time to molt into an
adult. At that time, the adult chews its way out of the cell. The length of time for pupal
development is dependent on caste (worker: 12 days, drone: 14.5 days, and queen: 8
days).

16

Figure 3: Life cycle of honeybees (Honeymell 2015). The life cycle of the honeybee begins with the
queen laying an egg which matures into larva, pupa, and finally an adult at day 21. The worker bees fulfill
their role by feeding the larva and sealing the cell.

Bee colonies display eusociality, an advanced level of organization characterized
by cooperative brood care, overlapping generations within a colony, and a division of
labor into reproductive and non-reproductive groups (Wilson 1971). With thousands of
honeybees inhabiting one hive, each class of bee performs a different job, a division of
labor described as temporal polyethism (Johnson 2008, Seeley 1985). The queen’s main
purpose is to mate and the drones’ is to inseminate the queen. The remaining jobs fall
onto the worker bees (Agriculture and Consumer Protection 1990). The nurse and forager
classifications of the worker bee are age-dependent. If the worker honeybee is three
weeks old or less, she is considered a nurse bee. Her roles are numerous and include:
cleaning the hive/comb, feeding the brood, caring for the queen, making orientation
flights, comb building and ventilating the hive, packing combs (with pollen, water, nectar
and honey), executions, and finally performing guard duty which is the last stage of the
nurse bee. The forager bee, a mature nurse bee, mainly concentrates on the needs of the
colony such as fetching nectar, pollen, water, and propolis (a resinous cement collected
by bees from buds of trees which is used as cement in repairing and maintaining the hive)
(Agriculture and Consumer Protection 1990, Kuropatnicki et al. 2013).
Due to their roles within the hive, the distinct classes of honeybees may vary in
their susceptibilities to environmental factors. A study performed by Vannette et al.
(2015) compared the expression of antimicrobial, immune and detoxification genes in
17

Apis mellifera between foragers and nurse bees. Researchers observed changes in
immune responses through ontogeny as the honeybee matures. In this specific example,
the forager bee exhibited greater expression of genes that were associated with immune
responses and detoxification activity than nurse bees. This heightened gene expression
was specifically pronounced in tissues that mediate nectar processing and social
interactions such as the mandibular gland and Malpighian tubules. The findings suggest
that the honeybee has mechanisms to deal with the different environmental threats that
could be encountered when out foraging. The ability of a specific gene’s expression to
shift as the honeybee makes the transition from nurse bee to forager is essential to its
function.
The Honeybee Genome:
Apis mellifera is the European honeybee, one of seven species and 44 subspecies
of honeybee (Engle 1999), that was used in this investigation on the effects of pesticides
on the honeybee transcriptome. The genus Apis is an ancient lineage of bees that evolved
in tropical Eurasia. The origin of Apis mellifera has been suggested to be Asia, the
Middle East or Africa. While the native range of Apis mellifera spans Europe, Africa, and
the Middle East, they are now found worldwide due to humans utilizing them as
pollinators and for their ability to make honey (Engle 1999). The karyotype of Apis
mellifera has been well characterized and consists of 16 chromosomes in haploid males
and 32 chromosomes in the diploid females (Figure 4). Apis mellifera lacks sex
chromosomes, a consequence of its haplodiploid sex determination.

18

Figure 4: Ideogram and karyotype for Apis mellifera (Honeybee Genome Sequencing Consortium
2006). The ideogram (in blue) shows the average chromosome lengths, positions, and sizes of the
heterochromatin bands. The percentage of heterochromatin mirrors the time of appearance of
heterochromatic bands. The karyotype is located below the ideogram.

The genome of Apis mellifera was published in Nature in 2006 by the Honeybee
Genome Sequencing Consortium, a partnership of over 250 scientists from 90 institutions
(Honeybee Sequencing Consortium 2006). This group established that honeybees have
approximately 11,000 genes, a surprisingly low number compared to other insects such as
Drosophila (13,600 genes) and Anopheles (14,000 genes) (Claudianos et al. 2006).
Possible explanations for this discrepancy are the bees’ haplodiploid system of sex
determination and their highly organized eusociality that limits the exposure of immature
bees to the external environment, lessening the need for functions that require
environmental interactions (Han et al. 2012).
a. Genome Organization
Honeybee genomes display several unique characteristics. (Honeybee Genome
Sequencing Consortium 2006). A honeybee’s genome can be distinguished from those of
other insects by possessing a high (A+T) content, high CpG content, and lacking major
transposon families. Due to the high A+T content, honeybee genes appear more
frequently in (A +T) rich domains. CpG is considered an over-represented dinucleotide
when compared to the mononucleotide frequencies. In genomes that have CpGs as the
target of cytosine methylases display a CpG deficit. The fact that the honeybee possesses
19

such high CpG content is surprising given that honeybees are the first protostomes to
possess DNA methylases, including CpG methyltranferase, that are similar to those found
in vertebrates. Methylated cytosines are known to be frequent sites of mutation, and the
expected Me-C

T mutations may contribute to the A+T richness observed in the

honeybee genome (Honeybee Genome Sequencing Consortium 2006). However, the
significance of cytosine methylation on the nucleotide composition of the honeybee
genome remains uncertain. The honeybee genome is also unique in that it lacks most of
the major families of transposons and retrotransposons. Lastly, the honeybee’s telomeres
appear to be relatively simple, lacking complex tandem repeats when compared to other
insects (Honeybee Genome Sequencing Consortium 2006).
b.

Functional Categories of Genes

The honeybee genome varies in both obvious and subtle ways compared to the
well-annotated Drosophila and human genomes. The honeybee genome shows greater
similarities to the vertebrate genome than the Drosophila and Anopheles genomes for
genes involving circadian rhythms, RNA interference, and DNA methylation. Honeybees
possess fewer genes than Drosophila and Anopheles for innate immunity, detoxification
enzymes, gustatory receptors, more genes for odorant receptors, and novel genes for
nectar and pollen utilization (Honeybee Genome Sequencing Consortium 2006). Key
examples of gene variation in honeybees in critical functional categories are described
below:
1. Immune System Pathways.
The eusocial lifestyle of the honeybee promotes favorable
conditions for infectious threats such as viruses, bacteria, fungi, protists,
and parasites. Honeybees protect themselves against such dangers by
grooming, which is considered a social defense, by raising young in
individual chambers within the hive, and by having the workers defend the
hive against potential vectors of disease (Honeybee Genome Sequencing
Consortium 2006). However, the honeybee genome, when compared to
other insect genomes, possesses fewer genes that are implicated in insect
immune pathways. While immune pathways such as Toll, Imd, and
JAK/STAT are intact, their functionality is reduced by two-thirds in
20

comparison to honeybee paralogues (Honeybee Genome Sequencing
Consortium 2006). This data suggests that honeybees utilize novel
immune pathways and possess immune systems that are focused on a
relatively small group of coevolved pathogens.
2. Anti-Xenobiotic Defense Mechanisms
Insecticide use is considered to be one of the factors that accounts
for the major loss of honeybee populations in parts of the world.
Detoxification enzymes are essential for honeybees to remove these
xenobiotics that they are exposed to in the environment. When compared
to Anopheles and Drosophila, the honeybee has approximately 30-50%
fewer genes that encode three superfamilies of xenobiotic detoxification
enzymes: carboxylesterase (CCE), cytochrome P450 (P450), and
glutathione-S-transferase (GST). These xenobiotic detoxification enzymes
are responsible for the metabolism of insecticides (Honeybee Genome
Sequencing Consortium 2006). Up- and down-regulation of these
different P450s can induce erratic behavior (Pereira et al. 2020). In other
insects, these three superfamilies have also been shown to be the frequent
source of mutations that confer insecticide resistance (Claudianos et al.
2006). The reduction of these anti-xenobiotic genes in bees may explain
their unusual sensitivity to insecticides.
A study performed by Mao et al. (2013) determined that pcoumaric acid, pinocembrin and pinobanksin 5-methyl ether (all
constituents found in honey) induce detoxification genes. By performing
RNA-seq analysis, it was shown that p-coumaric acid up-regulates all
classes of detoxification genes as well as some antimicrobial peptide
genes. This is of interest since p-coumaric acid could function as a
nutraceutical that regulates immune and detoxification processes. This
would explain why certain bees are not affected or less affected by toxins.
However, honey substitutes are often used which may compromise the
bee’s ability to deal with pesticides and pathogens thus contributing to
colony loss (Mao et al. 2013).
21

3. Antioxidants
Enzymes that break down free radicals are conserved in bees.
Unlike Drosophila and Anopheles, there is an expansion of the sigma class
of GST enzymes that protects against free radicals generated in aerobic
metabolism, and fewer numbers of the theta, delta, and omega classes that
protect from xenobiotics. Since free radicals are generated in bee flight
muscles during long foraging trips, selective pressure may be put on these
genes (Honeybee Genome Sequencing Consortium 2006).
4. Heat Shock Proteins
Heat shock proteins, or HSPs, are molecular chaperones produced
in cells in response to stressful conditions. Universal in all organisms, heat
shock protein genes are highly-conserved and assigned to families based
upon their sequence homology and molecular weight (Feder and Hoffman
1999). Honeybees have only 6 genes in comparison to the fly (5 genes)
and humans (8 genes). HSP70 proteins function in protein complexes with
HSP90 which are involved in signal transduction, ligand binding as well
as responding during cellular stress (Honeybee Genome Sequencing
Consortium 2006). Our laboratory adapted a HSP70 ELISA to monitor
cellular stress in different life stages of several species of bees (Barthell et
al. 2002, Hranitz and Barthell 2003, Hranitz et al. 2009a, Hranitz et al.
2009b, Hranitz et al. 2010).
Honeybees are one of few endothermic insects and have
astonishing levels of thermoregulation in colonies whose temperatures are
33-35°C. The core temperature of 33-35°C is maintained through bee
activities that warm or cool the colony as needed. While the hive
thermoregulates, the role of HSPs in the thermotolerance of individual
honeybees appears to function as in other animals, as honeybees do not
possess an increased number of genes that encode HSP70 family members
(Tong et al. 2019).

22

5. Circadian Rhythm and Behavior
Circadian rhythm of the honeybee is socially regulated (EbanRothschild and Bloch 2011). It is common to see differences in the
circadian clock among different castes. As a forager, the circadian clock is
used to anticipate day-night fluctuations in their environment, time visits
to flowers, dance language communication, and to determine sun-compass
orientation. The worker and the queen bee however have more flexibility
to their circadian rhythm. The internal clocks of the worker bees are
influenced by task specialization and regulated by the direct contact with
the brood. It is important to note that nurse bees, one of the youngest
bees in the caste system, do not possess circadian rhythms in behavior or
clock gene expression (Eban-Rothschild and Bloch 2011).
Sleep regulation is an important focus of the circadian clock.
Honeybees have a distinct sleep state with a characteristic posture,
reduced muscle tension, and elevated response threshold (EbanRothschild and Bloch 2011). The sleep state is considered a dynamic
process consisting of changes between deep and light sleep. A lack of
sleep often leads to an increase in the expression of sleep characteristics
the next day and will interfere with learning patterns (Eban-Rothschild
and Bloch 2011).
Period (per), timeless (tim), crytochrome (cry), clock (clk), cycle
(cyc), vrille (vri), and Par Domain Protein 1 (pdp1) are all considered
“clock genes.” The honeybee genome encodes a single orthologue for
each of these genes (Honeybee Genome Sequencing Consortium 2006).
The honeybee genome only encodes mammalian-type paralogues but does
not contain orthologues to Cry-d and Timeless1 genes that are considered
essential clock genes for Drosophila (Honeybee Genome Sequencing
Consortium 2006).
Sequencing of the honeybee genome revealed a cys-loop
neurotransmitter-gated ion channel superfamily. This superfamily contains
21 subunit members, two less than Drosophila, and an extra nicotinic
23

acetylcholine receptor subunit. The members within this superfamily are
known to contribute to honeybee behavior including foraging, learning,
memory, olfactory signal processing, mechanosensory antennal input, and
visual processing (Honeybee Genome Sequencing Consortium 2006).
Colony Collapse Disorder:
First described in France in 2006, Colony Collapse Disorder (CCD) was identified
as a potential cause for the rapid decline in honeybees and other pollinators around the
world. Honeybee populations in the United States have been steadily declining at a rate
of 1% per year since 1947 with steeper declines seen since 1987, with a majority of these
losses averaging 29-36% per year over the last four winters (Pilatic 2012). Figure 5
shows the United States managed honeybee colony loss within the past decade.

Figure 5: U.S. managed honeybee colony loss estimates (Bee Culture 2019). Acceptable levels
of winter loss are compared to total winter loss and total annual loss of honeybee colonies between
2006-2007 and 2018-2019. The total annual loss has surpassed the total winter loss and nearly
doubled, almost tripled, the acceptable levels for 2010-2011 through 2018-2019 years.

Symptoms of CCD are distinct from other loss epidemics and include (Pilatic
2012):
1. Colonies found suddenly empty of adult bees leaving their brood unattended.
2. No sign of dead bees.
3. No hive pests or food robbers despite there being plenty of honey and pollen
24

stores.
4. Common hive parasites are not present at levels that are thought to cause
population decline.
Most scientists agree that there is no one single cause of CCD, but rather a combination
of factors, such as nutritional stress, pathogens, and pesticides, that weaken bee colonies
by impairing the honeybee’s immunity (vanEngelsdorp et al. 2009).
Nutritional stress can affect colony health in several ways such as immune system
suppression and reducing reproductive viability of the honeybees (Pilatic 2012). A key
factor in nutritional stress is habitat loss. Due to the loss of habitat, the honeybee has a
less varied and nutritious diet (Pilatic 2012). Contributing factors may include a
decreased ratio of open to developed land and the use of broad-spectrum herbicides on
herbicide-resistant, genetically engineered crops (Benbrook 2009).
Pathogens have also been implicated as a cause of CCD. In many cases, this is
complicated by multiple pathogens present at the same time between separate bee
colonies (each colony expresses different diseases in various combinations) (du Rand et
al. 2015). The most common pathogens are parasitic mites, viruses and gut fungus.
Varroa mites, the most important pests to honeybees, rapidly increase in population size
after invading the hive and attach to developing larvae thereby devastating a colony. The
Varroa mite sucks hemolymph from both the brood and the adult bees, affecting bee
development within the brood and weakening the adult bees (Le Conte et al. 2010).
Varroa mites also act as vectors to transmit a number of viruses that can weaken the
colony (Zemene et al. 2015). Viruses that may contribute to CCD are deformed wing
virus and other related paralysis viruses that threatened bee survival by causing adult bees
to lose their ability to fly (Coulon et al. 2018). Also implicated in CCD is the fungus
Nosema, a pathogen that affects the digestive system of the workers (van Dooremalen et
al. 2018). Other potential reasons for CCD are new and emerging diseases, immunesuppressing stress due to any of the reasons previously mentioned, and a lack of genetic
diversity (United States Environmental Protection Agency 2016).
Pesticides as Environmental Stressors:
Pollinating insects are at risk for continuous exposure to a multitude of
agrochemicals (Pereira et al. 2020). In the U.S., over 1 billion pounds of pesticides are
25

used annually, amounting to $14 billion in sales (Atwood and Paisley-Jones 2017,
Alavania 2009). Over 1,200 active ingredients are approved, contributing to 18,000
different commercial formulations (Pilatic 2012). There are five major classes of
pesticides: organophosphates, carbamates, pyrethyroid, neonicotinyls and chlorinated
hydrocarbons (Agronomy 2017). A study performed by Frazier and colleagues (2008)
analyzed 108 pollen samples and identified a total of 46 different types of agrochemicals.
A total of 17 different types of agrochemicals were identified from a single pollen sample
in an Apis mellifera colony.
The neonicotinyls are a synthetic insecticide related to nicotine. Neonicotinoids
bind irreversibility to nicotinic acetylcholine receptors (nAchR), ligand-gated ion
channels on postsynaptic membranes. Acetylcholine is the major excitatory
neurotransmitter in the honeybee brain and controls a wide-range of behaviors essential
for survival. Due to the distinct features of their nAchR subtypes, insects are much more
sensitive to neonicotinoid pesticides than birds or mammals (Tomizawa et al. 2000).
Compared to other classes of insecticides, neonicotinoids are less toxic to birds and
mammals than insects leading to their widespread usage as insecticides (Agronomy
2017). Neonicotinoids have been an integral part of the global insecticide market since
the 1990s. At least 143 million of the 442 million acres of United States croplands are
planted with crops treated with 1 of 3 neonicotinoid pesticides that are known to be
highly toxic to bees (Pilatic 2012). These pesticides are clothianidin, imidacloprid and
thiamethoxam.
The mechanism of action of neonicotinoids on honeybee neurons has been
recently reviewed (Cabirol and Haase 2019). Neonicotinoids are acetylcholine agonists.
After binding to the nAchR, neonicotinoids produce a biphasic response. Initially, there is
an increase in the frequency of depolarization of the postsynaptic neuron followed by a
complete block to nerve signal propagation (Schroeder and Flattum 1984). Bees treated
with neonicotinoids initially show hyperactivity followed by convulsions and eventual
paralysis (Cabirol and Haase 2019). As shown in Figure 6, under normal conditions,
acetylcholinesterase in the synaptic cleft breaks down acetylcholine which prevents
overstimulation of the nAchRs. Neonicotinoids cannot be broken down by
acetylcholinesterase leading to overstimulation of the postsynaptic neuron. The specific
26

effects of neonicotinoids on the bee brain are complicated by the diversity of nAchRs and
their differential expression in various brain regions and developmental stages (Cabirol
and Haase 2019).

Figure 6: Neurotransmission with acetylcholine and neonicotinoids (Photo Credit: BioNinja NA).
Nicotinic acetylcholine receptors, also known as nAchRs, are cholinergic receptors. These receptors form
ligand-gated ion channels within the plasma membranes of postsynaptic neurons in the central nervous
system of the honeybee.

There are two classes of broad-spectrum neonicotinoids: nitroguanidines and
cyanoamines (Wu-Smart et al. 2016). Nitroguanidines are highly toxic to honeybees and
consist of imidacloprid, clothianidin, thiamethoxam and dinotefuran. Cyanoamines are
not as acutely toxic to the honeybee and include thiacloprid and acetamiprid (Wu-Smart
et al. 2016). Levels of these neonicotinoids in bee food sources are typically minute,
ranging from <1 to 8.6 ppb in nectar and <1 to 51 ppb in pollen (Goulson 2013)
(Appendix B-D). However, one type of exposure, corn seedling guttation, has been
shown to expose bees to neonicotinoid doses that exceed the LD50 (Girolami et al. 2009).
Even though neonicotinoid concentrations in nectar and pollen may be well below their
toxic lethal doses, they may produce sublethal effects in the bees that can have dramatic
effects on the colonies.
Environmental risks caused by neonicotinoids are still in question. Predominately
used as seed dressings, this process provides better targeting of the crop than spray
applications and with no action from the farmer, the crops are protected for several
months following sowing (Goulson 2013). Although there is evidence that neonicotinoid
application provides effective control of a range of insect pests, it is not clear as to the
27

extent this has on increased farm production or if there are increased economic benefits
when compared to alternatives (Goulson 2013). Other contributing factors to increase
production or profits are improved crop varieties, widespread use of artificial fertilizers,
new agronomic techniques and the development of successive generations of pesticides.
There is also a scant amount of studies to compare the effectiveness of neonicotinoids
with alternative means of pest control (Goulson 2013).
The persistence of neonicotinoids in soil is of great concern. Studies have shown
that between 1.6 and 20% of the active ingredient of seed dressings is absorbed by crops,
whereas the uptake by traditional spray applications exceed 50% (Goulson 2013). Up to
2% of the active ingredient is lost as dust during sowing, causing mortality to honeybees
flying nearby, and can be deposited on field vegetation at concentrations from 1-9 ppb
(Goulson 2013). More than 90% of the neonicotinoid active ingredient enters the soil and
can persist for 200, to > 1,000 days. Due to the high translocation rate of neonicotinoids
in plants, insecticides like imidacloprid, reach the flowers, leaving residue in the nectar
and pollen (Pereira et al. 2020). Foragers will continue to visit imidacloprid-treated crops
despite nectar and pollen contamination. These trips remain regular but at a slower rate
(Kessler et al. 2015, Pereira et al. 2020). Although the foragers’ cognitive process is not
significantly affected, the nectar odor that is brought back to the hive with them recruits
additional bees to visit the “contaminated” flowers, ultimately increasing the amount of
insecticide in the hive (Pereira et al. 2020).
The neonicotinoid imidacloprid became the world’s largest selling insecticide,
and second largest selling pesticide in 2008 with registered uses for over 140 crops in 120
countries (Jeschke et al. 2011). In 2010, it was estimated that 20,000 tons of imidacloprid
was produced globally, with 14,000 tons produced by China, which then exported 8,000
tons (Simon-Delso et al. 2015). With over 400 products for sale in the United States,
imidacloprid can be found as liquids, granules, dust, and packages that dissolve in water.
It is used for crops, inside and outside the home, as well as flea medications for animals
(Gervais et al. 2010). The two largest manufacturers of imidacloprid in the United States
are Crop Science and Albaugh LLC AgriStar (Figure 7; Bayer Crop Science 2017,
Albaugh LLC 2017).

28

C9H10ClN5O2

Figure 7: Chemical structure and molecular formula of imidacloprid (National Center for
Biotechnology Information 2020). IUPAC name of N-{1-[(6-Chloro-3-pyridyl)methyl]-4,5dihydroimidazol-2-yl}nitramide.

The number of deleterious effects that imidacloprid has on honeybees continues
to grow as more research is performed. Since neonicotinoids act on the central nervous
system, there are many physiological and behavioral impairments that have been
observed in worker bees. For example, bees exposed to sublethal doses of imidacloprid
have impaired foraging and homing abilities, suppressed immunity, delayed larval
development and reduced longevity, and diminished olfactory learning and memory
capacity (Henry et al. 2012, DiPrisco et al. 2013, Wu et al. 2012, Pereira et al. 2020,
Decourtye et al. 2004a, Decourtye et al. 2004b, Decourtye et al. 2003).
Pesticide Dynamics in Honeybees:
First defined by entomologist William Morton Wheeler in 1918, trophallaxis is
defined as the transfer of food or other fluids among members of a community
(Wainselboim et al. 2000). While worker bees acquire neonicotinoids directly from
consuming contaminated nectar or pollen, the queen bee may be exposed indirectly
through trophallaxis or food-sharing. A study by Wu-Smart and Spivak (2016)
demonstrated that environmentally relevant doses of imidacloprid had harmful effects on
queens (egg-laying and locomotor activity) and colony development (brood production
and pollen stores). These effects were less evident in larger colonies. Larger colony
populations may act as a buffer to pesticide exposure. Trophallaxis may help to lessen the
toxicity of pesticides by evenly distributing them among nest mates. This dilutes the
pesticides with uncontaminated food or bodily fluids already in the gut, rendering them
less harmful.
In addition to neonicotinoid pesticides, honeybees are exposed to other nonneonicotinoid pesticides such as fungicides and pyrethroids. A key mechanism used by
insects to counteract the effects of these toxins is metabolic resistance (du Rand et al.
29

2015). The major enzyme superfamilies that detoxify toxins are cytochrome P450
monooxygenases, glutathione transferases, and carboxylesterases. As previously noted,
there is a 30-50% or more reduction in the number of genes that encoded these enzyme
families compared to other insect genomes (Honeybee Genome Sequencing Consortium
2006). It is possible that because honeybees have fewer detoxification genes, their ability
to metabolize multiple toxins simultaneously are hindered (du Rand et al. 2015). For
example, fungicides inhibit cytochrome P450s that are important in insecticide
detoxification pathways (Johnson 2006). As previously mentioned, CCD has no single
cause which raises the possibility of synergy. For example, a pesticide in combination
with a fungicide can enhance toxicity. Studies have also shown synergistic effects
between pathogens and pesticides which led Poquet et al. (2015) to pose the question:
Does a pesticide promote pathogenicity of the pathogen or does the pathogen increase
toxicity of the pesticide?
The mechanisms used by healthy honeybees to detoxify neonicotinoid
insecticides are not well characterized. The goal of the du Rand study (2015) was to shed
light on the molecular mechanisms that honeybees use to detoxify nicotine, a toxic
natural alkaloid found in plants in the Solananceae family, such as tobacco. Nicotine and
synthetic neonicotinoid pesticides have similar modes of action; both mimic
acetylcholine by binding to nAchRs. This study used mass spectroscopy-based metabolic
and proteomic analysis to determine metabolic pathways and protein networks involved
in nicotine detoxification in newly emerged worker bees. The results indicated that
honeybees actively detoxify nicotine to its less-toxic metabolites, cotinine and cotinine
N-oxide (Phase I detoxification). This was followed by Phase II detoxification,
conjugation with glutathione catalyzed by glutathione-S-transferases. A total of 1470
proteins were identified with nearly 100 proteins that were up-regulated and 60 that were
down-regulated. The largest groups of up-regulated proteins include those whose
functions related to energy metabolism (ATP synthesis, glycolysis, TCA cycle enzymes,
etc.) suggesting increased energy production to support detoxification processes. Also,
up-regulated were proteins involved in detoxification, heat shock, and anti-oxidant
responses. The authors proposed that the increased energy production led to increases in
reactive oxygen species, which induced the expression of antioxidants and stress response
30

proteins. Based on their results, du Rand et al. (2015) proposed a model in which nicotine
exposure activates detoxification, oxidative, and stress response pathways in concert with
an increase in energy metabolism (Figure 8). These complex responses of honeybees to
nicotine provide a basis for us to predict the outcomes of our studies investigating the
effects of imidacloprid on the honeybee transcriptome.

Figure 8: Proposed mechanisms underlying the honeybee response to nicotine exposure (du Rand et
al. 2015).

Effects of Neonicotinoid Pesticides on the Honeybee Transcriptome:
While previous studies have found that sublethal doses of neonicotinoids impair
learning, memory capacity, foraging, and immunocompetence in honeybees, not much is
known about their molecular effects. Since undertaking this thesis work, several
laboratories have conducted studies to determine the effects of neonicotinoid pesticides
on the honeybee transcriptome. A study performed by Shi and colleagues (2017)
examined the transcriptome profile of honeybees after sub-chronic exposure to
31

thiamethoxam, a second-generation neonicotinoid, at 10 ppb over 10 days. Overall, there
were 609 differentially-expressed genes. Of these 609 genes, 225 genes were upregulated while 384 were down-regulated. The differentially-expressed genes were
mainly found to be associated with metabolism, biosynthesis, and translation (Shi et al.
2017). Zhiguo et al. (2019) reported down-regulation of brain genes related to immune,
detoxification, and chemosensory responses in honeybees after chronic oral exposure to
sublethal doses of imidacloprid. They proposed that this may contribute to decreased
olfactory learning abilities in imidacloprid-treated bees. Work by Wu et. al (2017) has
shown that genes encoding the major royal jelly proteins, essential for the health of
sustainable colonies, were strongly down-regulated in honeybee larvae exposed to
sublethal doses of imidacloprid.
Effects of neonicotinoids on the honeybee brain appear to be dose dependent.
Christen et al. (2018) reported that honeybees treated for 48 hours with 0.3 ng/bee of
imidacloprid had 7 coding regions up-regulated and 19 regions down-regulated. When
this study was repeated with 3 ng/bee of imidacloprid, a total of 36 regions were upregulated and 77 regions were down-regulated providing evidence that different genes are
transcribed when bees are exposed to different amounts of imidacloprid. In particular,
genes related to metabolism and detoxification were differentially expressed, generally in
a concentration-dependent manner. Collectively, these studies highlight the utility of
transcriptome profiling in providing insights into the mechanisms mediating the toxicity
of pesticides. Clearly, future studies are needed to tease out such details as dose
dependency, chronic versus acute exposures, laboratory versus field exposures,
synergistic effects, and differences among caste members and developmental stages.
Importantly, for a better understanding of the effects of neonicotinoids on honeybees,
future studies are needed that link their molecular effects to physiology and behavior.
Sublethal Stress in Honeybees:
Regulatory agencies assess the impact a pesticide or chemical agent has by
determining the acceptable risk. Acute toxicity is quantified by determining the dose at
which half of the insects die within a specific time frame, denoted as LD50 as seen in
Figure 9 (Pilatic 2012).

32

Figure 9: Dose-response relationship for acute toxicity (Möller NA). Also known as lethal toxicity, the
LD50 (denoted by the grey arrow) is determined by the dose at which half of the insects die within a specific
time frame.

The symptoms of acute toxicity include agitation, vomiting, wing paralysis, arching of
the abdomen, and uncoordinated movement. More common in field exposure, sublethal
toxicity has symptoms of disorientation, reduced mobility and foraging, impaired
memory and learning, and shifts in communication behavior (Blacquière et al. 2012).
While these sublethal effects do not kill individual bees, they may have a profound
impact on the dynamics and function of the whole colony. For example, these sublethal
impacts may interfere with the food collection process for the colony impairing foraging,
the proboscis extension reflex, olfactory learning and memory (Blacquière et al. 2012).
Sublethal effects may also negatively impact reproductive rate, survival of larvae, and
recruitment into nurse castes which ultimately weakens the hive, leading to colony losses
that may not become apparent for weeks (Blacquière et al. 2012). Furthermore, the
reduced viability and quantity of the stored sperm in the mated queen may compromise
the queen’s breeding success (Wu-Smart and Spivak 2016).
A leading hypothesis for CCD is sublethal stress due to environmental stressors
such as pesticides, disease and parasites, habitat change and loss (Bryden et al. 2013).
Redundancy in social bee colonies allows for a significant loss of their workers without
any significant impact on colony function and productivity. However, if bees become
33

impaired rather than die, this may impose undue stress on the colony which can lead to a
cumulative effect on normal colony function through indirect mortality and production
losses not attributed to pesticides (Bryden et al. 2013). The Sublethal Stress (SLS) Model
introduced by Bryden et al. (2013) incorporates the effects of sublethal stress on colony
function to explain these colony dynamics. The SLS Model (Figure 10) predicts that
colonies can either persist or become extinct, dependent on the initial conditions that
determine whether or not the colony exceeds critical reproductive and mortality rates
(Bryden et al. 2013).

Figure 10: Sublethal Stress (SLS) Model (based on the Bryden model 2013). The model shows the
relationship between the intensity of a stressor and its duration and how this could lead to no effect,
sublethal effects or mortality.

To test this mathematical model, Bryden’s laboratory compared colony dynamics
in bumble bees treated with field-relevant doses of imidacloprid. The patterns of colony
sizes, birth rates, and death rates fit the SLS Model. Bumble bee colonies failed when
they were exposed to sublethal levels of imidacloprid for an extended period of time
resulting in a decrease in colony function (Bryden et al. 2013). By testing models against
data collected from failing colonies, it was concluded that social bee colonies have
positive density dependence, are subject to an Allee effect (that population size or density
is correlated with the mean individual fitness of a population), and that there is a critical
stress level for the success of a colony (Bryden et al. 2013, Drake and Kramer 2011).

34

They concluded that small increases in the level of stress can swing the pendulum
towards success or failure of the hive.
The SLS Model was compared to two alternative models: the Khoury Model and
the LA Model. The Khoury Model included lethal stress but not the impairment and
feedback caused by sublethal stress whereas the LA Model includes the toxic effects
from pesticides at the larval stage (Bryden et al. 2013). These three models were fitted
using the NISS algorithm that calculates a likelihood value for a model based on all
possible paths. From this modeling, it was determined that the SLS Model best matched
the pattern of birth rates decreasing and death rates increasing in the treated colonies. It
was concluded that colony function is important, in explaining the dynamics of the
treated colonies and therefore suggests a mechanism by which sublethal effects on
individual bees can lead to colony failure (Bryden et al. 2013). This study demonstrates
that sublethal stressors must have a chronic impact before effects are observed and that
stressors that impair colony function cause an Allee effect, making colonies susceptible
to stress at earlier points in their life cycle. The authors point out the irony that the
elaborate social organization that leads to the success of social bees may also be a key
factor contributing to colony failure and their population declines.
Cellular Responses to Stressors by Honeybees:
A general indicator of cellular stress is Heat Shock Protein 70 (HSP70). Although
normally present in cells, HSP70 levels become elevated at times of stress to help
maintain stress resistance (e.g., Hranitz et al. 2009). HSP70, a molecular chaperone,
works by binding peptides to prevent misfolding during exposure to foreign toxic
substances or other stressful conditions. Once the stressor is removed, the peptides are
released and are able to return to their normal cell functions. If stress levels become too
elevated, it is possible for HSP70 to activate apoptotic mechanisms and cause cell death
to avoid an inflammatory response. When cells are exposed long-term to stress,
irreversible consequences such as loss of nervous system control, delayed growth, and
inability to process sensory input for foraging may occur (Feder et al. 1997, Jones et al.
2007). HSP70 is an excellent marker for measuring cellular stress in honeybees and has
proven useful in a variety of applications, such as evaluating management practices,
assessing seasonal changes, and in toxicological studies (Hranitz et al. 2009).
35

Although there are linear and threshold models of stress, the hormesis toxicology
model can explain the entire range of physiological stress responses. Hormesis is defined
as a dose-response relationship that is stimulatory at low doses and inhibitory at higher
doses (Deng et al. 2001). The hormetic model illustrates variable stimuli proportions as
parabolic curves. This allows for inferences on positive, negative, and peak response
levels across a range of treatments (Calabrese 2008). Figure 11 illustrates hormesis as it
relates to stress and effect or dose and response.

Figure 11: Hormesis dose-response curve (Merritt 2011-2020). Hormesis is a U- or J-shaped doseresponse curve characterized by a low-dose stimulatory and high-dose inhibitory responses. A stimulus that
produces a harmful biological effect at a moderate to high doses may produce beneficial effects at lower
doses (Calabrese 2008).

The hormesis model was useful in describing the response of honeybees to
ethanol in a study performed by Hranitz and colleagues (2010). This study, based on
work by the Abramson Laboratory at Oklahoma State University, previously
demonstrated that ethanol intoxication has dramatic effects on learning and behavior in
the bee. Intoxicated honeybees were less active, had poor motor coordination,
demonstrated preference for sugar in ethanol solutions, showed increased levels of
aggression, displayed impaired foraging decisions and poor communication, and had a
similar time course of elevated blood alcohol elevation as humans (Abramson et al. 2000,
2002, 2003, 2004a, 2004b, 2005; Bozic et al. 2006, 2007). Levels of HSP70 were
measured at 4-hours post-ingestion in control bees fed 1.5 M sucrose or bees fed with
2.5%, 5%, or 10% ethanol (Figure 12; Hranitz et al. 2010). A hormetic response was
36

observed with peak levels of HSP70 noted at 5% ethanol. HSP levels were lower at
higher (10%) and lower (2.5%) ethanol concentrations.

Figure 12: HSP70 concentrations among pretreatment groups of honeybees and those treated with
ethanol (Hranitz et al. 2010). The above graph shows high stress levels, as indicated by HSP70, appeared
when cellular stress was induced with ethanol intoxication. A hormesis curve is noted among the treatment
group. A is significantly different than B at p < 0.05.

RNA microarray studies of gene expression in the honeybee brain following the fourhour exposure to ethanol showed significant changes in 609 genes. Gene ontology
analysis revealed changes in brain metabolism, cellular stress, and signaling pathways,
protein synthesis, and carbohydrate and lipid metabolism (Hranitz et al. unpublished). A
similar approach will be applied to this investigation to determine the effects of sublethal
doses of imidacloprid on the honeybee brain transcriptome.
Oxidative Stress:
Pesticides have been found to produce oxidative stress in honeybees (Henry et al.
2005). Oxidative stress occurs when reactive oxygen species accumulate in an organism
and cause damage (Berlett and Stadtman 1997). Reactive oxygen species (ROS) are
unstable chemical species containing oxygen that are formed during aerobic metabolism.
Examples of ROS include peroxide, hydroxyl radicals, and superoxide. Under normal
circumstances, ROS are in low concentrations. When an organism is under high stress,
ROS begins to accumulate in cells causing damage to DNA, lipids, and proteins. ROS
have been linked to diseases such as cancer, diabetes, and heart disease in humans. ROS
37

are degraded by antioxidant enzymes such as superoxide dismutase, catalase, and
glutathione transferase. Superoxide dismutase (SOD) converts superoxide, a powerful
ROS, to the less dangerous hydrogen peroxide. Catalase and glutathione transferase
(GST) work to convert hydrogen peroxide to water. High levels of these antioxidant
enzymes in organisms is an indicator of oxidative stress (Finkel and Holbrook 2000).
The link between pesticides and levels of antioxidants enzymes in honeybees has
been demonstrated in recent field studies. A study performed by Dussaubat et al. (2016)
reported increased catalase and GST activity in the head of queen bees following
exposure to environmentally relevant concentrations of imidacloprid. Chakrabarti et al.
(2014) compared the effects of the pesticide paraquot on the levels of antioxidant
enzymes in laboratory and field populations of two native Indian bees. The results
showed elevated levels of antioxidant enzymes in bees exposed to sublethal doses of
pesticides in the field and in the laboratory compared to controls. The authors postulated
that these higher levels of antioxidants help to protect bees during exposure to
environmental stressors. Oxidative stress has been shown to decrease the survival and
homing abilities of honeybees (Simone-Finstrom et al. 2016). This is important since
most crops are pollinated by honeybees that are transported from state to state, a
procedure called migratory management. It has been found that migratory management
induces oxidative stress in bees and produces effects that are similar to that of CDD. The
study concluded that there was a significant decrease in lifespan and foraging capabilities
of the affected bees (Simone-Finstrom et al. 2016). Together these studies underscore the
need to understand the mechanisms underlying the toxicity of neonicotinoid pesticides
and their contribution to pollinator decline.

II.

Objectives and Research Questions

Goal: The overall goal of our research is to describe the integrated responses by the
honeybee to sublethal doses of the common neonicotinoid, imidacloprid. This research
expands on previous studies in our laboratory, and on independent studies, by
investigating the physical, cellular, and molecular responses to imidacloprid by the
honeybee brain.

38

Specific Objectives: This research was conducted in two phases- Experiment 1 and
Experiment 2. The specific objectives of each are:
•

Experiment 1: The Effects of Sublethal Doses of the Neonicotinoid
Pesticide Imidacloprid on Motor Function and Cellular Responses in
Honeybees.
o To determine the effects of sublethal doses of imidacloprid on
motor function and overall cellular stress as determined by the
levels of Heat Shock Protein 70 and the oxidative stress enzyme
Superoxide Dismutase.
o To identify a sublethal dosage of imidacloprid that yields a
significant cellular stress response that will be used in Experiment
2.

•

Experiment 2: The Effects of Sublethal Doses of Imidacloprid on Gene
Expression in the Honeybee Brain
o To compare overall gene expression in brain tissue of the control
and imidacloprid treatment groups.
o To determine which functional classes of genes are up-regulated or
down-regulated following imidacloprid treatment.
o To examine differences in gene expression that occur following
imidacloprid treatment in pathways regulating detoxification, heat
shock proteins, oxidative enzymes, energy metabolism, circadian
rhythms, cell signaling, apoptosis, and longevity.

Hypotheses:
•

Bees exposed to sublethal doses of the neonicotinoid imidacloprid will display
impaired motor functions and a significant cellular stress response reflected by
increased levels of Heat Shock Protein 70 and the oxidative enzyme Superoxide
Dismutase.

•

A dose of imidacloprid can be determined that yields a significant cellular stress
response in honeybee brains.

•

Gene expression patterns will be altered in the brain tissue in imidacloprid treated
bees compared to controls.
39

•

Bees treated with imidacloprid will express different functional classes of genes
compared to control bees.

•

Pathways regulating detoxification, heat shock proteins, oxidative enzymes,
energy metabolism, and apoptosis will be up-regulated, while pathways regulating
circadian rhythms, cell signaling, and longevity will be down-regulated in
imidacloprid-treated bees.

Significance of the Study: This study investigates how acute sublethal exposure to
imidacloprid, the most widely used neonicotinoid pesticide, affects gene expression in an
economically crucial pollinator. An understanding of the specific gene networks and
cellular pathways affected by sublethal imidacloprid intoxication may help the scientific
community to better understand the mechanisms of neonicotinoid toxicity underlying the
sublethal pesticide effects contributing to pollinator declines. This work can serve as a
springboard for future hypothesis-driven gene expression studies that relate specific
molecular changes to biological functions and organism-level performance. This
integrated approach, connecting the responses of organisms in the field to their
underlying cellular and molecular mechanisms, is essential to understanding and
preventing CCD.

III.

Methods

Experiment 1: The Effects of Sublethal Doses of the Neonicotinoid Pesticide Imidacloprid
on Motor Function and Cellular Responses in Honeybees
a. Experimental Design
Figure 13 outlines the overall design of Experiment 1. First, honeybees were
collected, harnessed, and fed with 1.5 M sucrose to satiation. Surviving, healthy
honeybees were randomly assigned to control or treatment groups after 22-24 hours after
harnessing and feeding. The total number of honeybees is n=149 (Negative Control= 20,
Positive Control= 18, 1/5th= 19, 1/10th= 18, 1/20th= 17, 1/50th= 20, 1/100th= 20, and
1/500th= 17). Control groups were fed 1.5 M sucrose while treatment groups received
sublethal doses of imidacloprid in 1.5 M sucrose (Appendix E and F). After 4 hours, bees
were tested for motor responses and bee heads were removed and frozen for subsequent
40

measurement of HSP70, a marker of cellular stress, and SOD, a marker for oxidative
stress. After evaluating these results, an optimal sublethal dose of imidacloprid was
selected to dose bees in Experiment 2 for determining its effects on the bee brain
transcriptome.

Figure 13: Experiment 1 design. Upon collection, the honeybees were harnessed in modified 1.5 mL
microcentrifuge tubes and fed to satiation. After 22-24 hours at room temperature (22°C), bees were
randomly assigned to a control or treatment group. At 4 hours after treatment, motor scores were assessed,
or bee heads were removed for analysis of HSP70 and SOD.

b. Honeybee Collection:
Honeybees were collected locally at hives maintained by Dr. John M. Hranitz,
Bloomsburg, PA. Bees were recruited to a feeder located 10 meters from the hive. The
feeder consisted of an inverted mason jar containing a 50% sucrose solution with a few
drops of orange extract (Figure 14a). Bees were then captured in groups of 2-3 in clear
plastic collection jars with holes in the lid and transported to the laboratory on campus in
a cooler (Figure 14b). Bees were anesthetized at -20°C for about 3-5 minutes until
immobile and then were restrained in harnesses modified from 1.5 mL microcentrifuge
tubes with thin strips of duct tape placed between the head and thorax (Figure 14c). Bees
were fed 1.5 M sucrose to satiation and held for 22-24 hours at room temperature (22°C).

41

Figure 14: Honeybee feeder, collection, and harnessing. A) A glass mason jar inverted on the lid of a
petri dish was used as a feeder. The feeder contained a 50% sucrose solution and orange extract as an
attractant. B) Once a significant number of bees were attracted to the feeder, carefully, two or three bees
were captured into a clear collection jar with a lid that contained holes so that they could be transported. C)
Honeybees were harnessed in modified 1.5 mL microcentrifuge tubes and secured with a thin piece of duct
tape.

Negative and positive control bees were fed 1.5 M sucrose. The bees in the
negative control group remained at room temperature (22-25°C) while bees in the
positive control groups were placed in an incubator set to (42°C). A temperature of 42°C
was chosen because it falls slightly below the critical thermal maximum for Apis
mellifera, 42.8 ± 2.8°C (Atmowidjojo et al. 1997). Cellular stress responses have been
shown to reach their peak slightly below the critical thermal maximum (Atmowidjojo et
al. 1997). Honeybees in the treatment groups were fed with imidacloprid (Macho®4.0,
AgriStar, 40.07 g imidacloprid/100 g solution) in 1.5 M sucrose at doses of 3.6 ng/μL
(1/5th), 1.8 ng/μL (1/10th), 0.9 ng/μL (1/20th), 0.36 ng/μL (1/50th), 0.18 ng/μL (1/100th)
and 0.036 ng/μL (1/500th) of the LD50 (18 ng/bee) (Figure 15; Karahan et al. 2015). Both
control and treatment groups were monitored for four hours. At the completion of the
treatment, motor coordination was evaluated, and bee head capsules were removed and
frozen at -80°C to be used in HSP70 or SOD assays.

42

Figure 15: Honeybee treatment. Honeybees were randomly assigned following the satiation period. The
treatment groups were color coordinated, labeled with the treatment doses, and numbered. The honeybees
were separated and faced away from each other to avoid trophallaxis.

c. Motor Testing:
Four hours after treatment, all groups of honeybees (n=149) were scored in four
categories of motor function: leg movement, abdomen movement, antennae
responsiveness, and proboscis extension reflex. Leg and abdomen movement were simply
observed while the bee was in the harness. The antennae and proboscis extension reflex
were tested by bringing a Q-tip dipped in 1.5 M sucrose close to the antennae and
proboscis and observing the response. Each test was scored either 0 for no function, 1 for
impaired function, or 2 for normal function. The maximum score that each bee could
obtain was 8 and a minimum score was 0.
d. Homogenization of Bee Head Capsules
After treatment, frozen head capsules were homogenized in microcentrifuge tubes
with a pestle in phosphate buffered saline (PBS) containing 0.2% sodium azide, 2 mM ptosyl-L-arginine methyl ester (TAME), and a protease inhibitor cocktail (Roche), pH 7.6.
The homogenates were then centrifuged at 16,000 g for 20 minutes at 4°C. The
supernatant was removed and stored at -80°C.
e. Superoxide Dismutase (SOD) Assay
Superoxide Dismutase Colorimetric Activity Kit by Arbor Assays (Kit K028-H1)
was used to test for SOD in thawed bee head supernatants. SOD neutralizes superoxide
radicals (O2-). In this assay, xanthine oxidase produces superoxide in the presence of
oxygen. The superoxide produced converts a colorless substrate into a yellow-colored
product. The reactions involved in the assay are shown below (Figure 16). When more
43

SOD is present in the samples, the superoxide concentration will decrease, resulting in
less yellow-colored product.

Figure 16: Superoxide dismutase reaction and assay principle (Arbor Assays).

Standards were prepared by making serial dilutions of the bovine erythrocyte
SOD standard provided in the kit at concentrations of 4.0 U/mL, 2.0 U/mL, 1.0 U/mL,
0.5 U/mL, 0.25 U/mL, 0.125 U/mL, and 0.0625 U/mL. Ten µL of the appropriate
standards and samples were pipetted in duplicate into a 96-well Corning Costar 3695
plate. Ten µL of assay buffer was pipetted into the top row of wells to serve as a blank.
This was followed by additions of 50 µL of substrate solution and 25 µL of xanthine
oxidase solution to each well. The plate was incubated at room temperature for 20
minutes and the absorbance was then read at 450 nm using a Tecan Genios Plus
Microplate Reader and Magellan software (Figure 17).

Figure 17: SOD assay microplate. Xanthine oxidase produces superoxide in the presence of oxygen. The
superoxide produced converted the colorless substrate into a yellow-colored product. Greater levels of SOD
in the sample result in less of a yellow-colored product.

44

SOD activity was calculated for each sample in U/mL using MyAssay software. Samples
were corrected for protein content by dividing the SOD activity by the protein
concentration in the sample which had previously been determined.
f. Protein Assay
The protein concentration of each sample was determined using a colorimetric
dye-binding assay (Bio-Rad RC DC Protein Assay Kit). The homogenization buffer
served as blanks and standards were prepared by making serial dilutions of Bovine Serum
Albumin (BSA stock- 7.1 μg/μL) in homogenization buffer (Appendix G).
Using the Bio-Rad Protein Assay Kit, Reagent A1 was prepared by adding 50 μL
of Reagent S to 2.5 mL of Reagent A. Five μL aliquots of the protein standards, blanks,
and samples were pipetted in triplicate into a 96-well Corning Costar 3695 plate. Next, 25
μL of A1 was added to each well, followed by 200 μL of Reagent B. Plates were then
gently agitated using a rotator. After 15 minutes, absorbance was measured at 750 nm
using the Tecan Genios Plus Microplate Reader (Figure 18).

Figure 18: Bio-Rad DC Protein Assay. The variations in the intensity of the blue color represent different
amounts of protein in the samples. Wells remained a greenish-yellow color if there was no protein or very
little protein in the sample.

g. HSP70 ELISA
HSP70 levels were quantified using monoclonal ELISA with primary antibody
Mouse anti-Bovine HSP70 (Sigma H5147) and secondary antibody Goat anti-Mouse IgG
Horseradish Peroxidase (Sigma A0168) (Barthell et al. 2002, Hranitz and Barthell 2003).
It has been shown by immunoblotting that the primary antibody binds both the
constitutive or cognate and the inducible forms of HSP70 (Sigma).

45

Using the results obtained from the BioRad Protein Assay, thawed homogenate
supernatants were diluted to a protein concentration of 400 ng/μL in homogenization
buffer and stored at -80°C. A standard curve was prepared using stock bovine HSP70
(1μg/μL) (Sigma H9776) to yield final amounts ranging from 10 to 500 ng per well.
Homogenate buffer served as blanks. Five μL of homogenization buffer, HSP standards
(Appendix H), or samples were loaded to designated cells in triplicate. Each sample well
contained equal amounts of soluble protein, 2000 ng. 195μL of binding buffer (10 mM
Na2CO3/NaHCO3, pH 9.6) was added to each well using a multichannel pipettor, covered
in plastic wrap, and incubated at 4°C overnight.
After 24 hours, the binding buffer was removed from the microplate by inversion.
Each well was washed once with 200 μL of Phosphate Buffered Saline containing 0.05%
Tween 20, pH 7.6 (PBST). The PBST was removed by inversion. Known as the blocking
step, 200 μL of PBST with 1% BSA was added per well, wrapped in plastic wrap, and
incubated at 37°C for 1 hour. The wells were then washed once with 200 μL of PBST.
Mouse anti-bovine HSP70 (primary antibody) was added (200 μL per well) and
incubated at 37°C for one hour. Plates were subsequently washed four times with 200 μL
of PBST. The secondary antibody (Goat Anti-mouse Horseradish Peroxidase) was added
(200 μL per well) and incubated at 37°C for one hour. Each microplate was then washed
six times with 200 μL of PBST. Following the series of washes, 150 μL of TMB
Substrate (BioRad) was added per well, the plate covered with plastic wrap, and
incubated at room temperature on a shaker for up to 30 minutes. During this time, the
samples turn various shades of blue. To stop the reaction, 100 μL of 1M sulfuric acid was
added resulting in different shades of yellow (Figure 19).

Figure 19: Monoclonal antibody ELISA test for HSP70. The amount of HSP70 in each well is
determined by the intensity of the yellow color produced when sulfuric acid was added to halt the reaction.

46

The absorbance was then read at 450 nm using the Tecan Genios Plus microplate reader.
The slope of the HSP Standard Curve was used to calculate the HSP concentration of
each sample.
h. Statistics:
Motor function tests were analyzed with a nonlinear regression fitted using the
sigmoidal dose response curve in GraphPad Prism v. 6.0 was performed. Chi square
analysis was performed on the frequencies (counts) of bees in each category for each
treatment (SPSS Statistics). Values of p<0.05 were considered significant. For SOD and
HSP70 Assays, one-way ANOVA and Post-Hoc Tukey tests were conducted using JMP
Pro v.12 software to test for difference among the treatments. To construct the TD50
curve, the percentage of impaired bees (motor score <8) was graphed versus the dose of
imidacloprid received. Values of p<0.05 were considered significant. There was no
correction for multiple testing for motor function, SOD, and HSP70 performed.

Experiment 2: The Effect of Sublethal Doses of Imidacloprid on Gene Expression in the
Brain
a. Experimental Design
The overall design of Experiment 2 is shown in Figure 20. Honeybees were
collected, as described in Experiment 1, and were randomly assigned to four groups:
Control-0h (C0), Control-4h (C4), Treatment-0h (T0), and Treatment-4h (T4). Control
groups were fed 10 μL of 1.5 M sucrose. Treatment groups received 1/20th of the LD50
(0.9 ng/bee) of imidacloprid in 10 μL of 1.5 M sucrose. This sublethal dose was selected
based on the results of Experiment 1 as one of the peak values in the hormetic HSP70
response to the pesticide but one dose lower than the dose that produced impaired
locomotor function. Honeybee heads were excised at time 0 or 4 hours and the brains
dissected. Total RNA was isolated from the brain tissue and transported overnight on dry
ice to the Keck Center at University of Illinois Urbana-Champaign. Genome-wide RNA
sequencing was conducted by the Keck Center. Following post-sequencing analysis,
pathway analysis was used to identify pathways altered by imidacloprid treatment.

47

Honeybee
Collection

Experiment
2
Treatments

Isolate Total
RNA

RNA
Sequencing

PostSequencing
Analysis

Pathway
Analysis

Figure 20: Experiment 2 design. Roughly 10-14 honeybees were assigned to each group and treated
(Control-0h, Control-4h, Treatment-0h, and Treatment-4h). Bee heads were removed, and total RNA was
isolated in dissected brains. The isolated RNA was sent to the University of Illinois for RNA sequencing.
Lastly, pathway analysis was performed to identify key pathways of interest.

b. Honeybee Collection and Treatments:
Honeybees were collected using the same method as described in Experiment 1.
Bees were harnessed and fed 1.5 M sucrose to satiation. After 22-24 hours, the bees were
randomly assigned to Control-0h (C0), Control-4h (C4), Treatment-0h (T0), and
Treatment-4h (T4). Treatment groups received 1/20th of the LD50 (0.9 ng/bee) of
imidacloprid that was determined based on the results of Experiment 1. A 1/20th LD50
dose of imidacloprid was a sublethal dose that was enough to elicit a stress response
without being debilitating. Upon completion of treatment (0 Hour or 4 Hour), the bee
head capsules were removed using sterile scissors and stored at -80oC in a sterile tube
until brain dissection was performed. The total number of honeybees collected was n= 43
(C0= 12, C4= 11, T0= 9, T4= 11).
c. Head Capsule Dissection:
Head capsules were removed from the -80oC freezer and placed on dry ice. The
dissecting microscope, microscalpel, forceps, petri dish, and surrounding lab bench was
cleaned with RNase Away (Molecular BioProducts). The honeybee head was placed in a
petri dish containing dry ice and 200 proof ethanol (Figure 21a). The anterior exoskeleton
was removed with a microscalpel exposing the brain (Figure 21b). Further dissection
removed retinal tissue and the compound eye lens (Figure 21c). Finally, the bee brain was
removed from the rest of the remaining head capsule (Figure 21d).

48

Figure 21: Procedure for isolating the honeybee brain (Photo Credit: Megan Duell 2014). A) Bee head
capsules were placed on a slurry of dry ice and ethanol. B) The exoskeleton is chipped away during
dissection, exposing the brain. C) Eye tissue and the lens are removed to see the boundary between retinal
and brain tissue. D)The ivory-colored brain resembles the shape of a bat and is removed. The blue arrow
indicates the actual size of the excised brain shown on the tip of a small pestle in comparison with fingers
in the background.

d. RNA Isolation
RNA was isolated from honeybee brains using a modified Trizol-RNeasy method
based on the procedure outlined in the Qiagen RNeasy kit as described by Sen Sarma et
al. (2009). The Qiagen RNeasy technology for purifying RNA, outlined in Figure 22,
combines the selective binding of RNA to a silica-based membrane with a fast microspin
technology.

49

Figure 22: Schematic of RNA isolation procedure (Qiagen 2020). RNeasy Protect Mini Procedure for
animal tissue is illustrated above. The procedure was modified to perform this research.

Standard procedures for handling RNA were employed during the RNA isolation
steps. All surfaces and instruments were cleaned with RNase-Away (Molecular
BioProducts) to remove RNases that are abundant in the environment. Gloves were worn
at all times. After homogenization with TRIzol (Invitrogen), the specimens were
incubated on the benchtop for five minutes. Equal parts of RNase free water and
chloroform were added, vortexed for 15 seconds, and incubated for an additional three
minutes. The specimens were then centrifuged at 14,500 RPM for 15 minutes at 4°C.
Following centrifugation, the aqueous phase was carefully removed without
disturbing the interface and placed into a sterile RNase Free 1.5 mL collection tube. An
equal amount of 70% ethanol was added to each sample and vortexed for 10 seconds.
Approximately 700 μL of the aqueous-ethanol solution was transferred into a RNeasy
50

mini spin column sitting in a 2 mL collection tube and centrifuged for 30 seconds at
12,000g at room temperature. The flow-through was poured back into the column which
was placed into a new 1.5 mL collection tube and centrifuged again for 30 seconds at
12,000g at room temperature discarding the flow-through after centrifugation was
complete.
Forty μL of a DNase/RDD buffer stock mixture was added to the column and
incubated at room temperature for 15 minutes. After this, 350 μL of RW1 buffer was
added to the RNeasy column and centrifuged at 12,000g for 15 seconds at room
temperature. Once the column was transferred to a new 2 mL collection tube, 500 μL
RPE buffer was added onto the column and centrifuged for 15 seconds at 12,000g at
room temperature. The liquid was decanted from the tube and the column was replaced.
An additional 500 μL RPE buffer was added to the column and centrifuged for 2 minutes
at 12,000g at room temperature. The liquid was decanted, and the column replaced. The
samples were then centrifuged at 12,000g for one minute to remove any residual ethanol.
The RNeasy column was transferred to a new 1.5 mL tube. Thirty μL of RNase free
water was pipetted onto the column and centrifuged at 12,000g for one minute. An
additional 30 μL of RNase free water was added onto the column and centrifuged for one
minute. The purified RNA in the eluate was stored in 1.5 mL microcentrifuge tubes at 80°C.
The quality and quantity of the purified RNA was assessed by determining the
ratio of absorption at 260 nm and 280 nm using the Nanodrop (Spectrophotometer ND1000) to confirm the amount and quality of the RNA present. For high-quality RNA, the
desired A260/A280 ratio is in the range of 1.9-2.1 (Mater Methods 2012). Table 1
compiles the average absorption ratio and the average RNA concentration along with the
calculated standard deviation for each testing group.

51

Table 1: Average absorption ratio and average protein present. The average A260/280 ratio and the
average RNA concentrations were determined for each treatment group (n=6).

Treatment
Group

Average
Concentration
A260/A280
Ratio (ng/μL)

Control-0h
Control-4h
Treatment-0h
Treatment-4h

1.88
1.69
1.90
1.86

Concentration
A260/A280
Ratio (ng/μL)
Standard
Deviation
0.069
0.061
0.099
0.043

Total RNA
(ng)
Average

Total RNA
(ng)
Standard
Deviation

38.0
42.3
48.2
32.8

3.92
11.9
9.05
10.9

The six highest quality RNA samples for each treatment group were selected based on the
A260/A280 ratio for RNA sequencing. To confirm RNA presence, a formaldehydeagarose gel electrophoresis was performed on four specimens chosen at random.
Formaldehyde was included in the gel to denature RNA and inhibit RNAses (Figure 23).

Figure 23: Formaldehyde gel for RNA. Formaldehyde-agarose gel electrophoresis was run to confirm the
presence of RNA. Ribosomal RNA bands are denoted by the blue arrows and mRNA is a range of sizes
indicated as a smear. The formaldehyde gel contained: 1% agarose gel, MOPS (20 mM MOPS(3-(Nmorpholino) propanesulfonic acid), 2 mM sodium acetate, 1 mM EDTA, pH 7.0), and 2.2 M formaldehyde.
Loading dye contained 2 microliters of ethidium bromide. The gel was run for approximately 45 minutes at
95 volts.

The purified RNA samples were shipped to the Keck Center for Comparative and
Functional Genomics at the University of Illinois Urbana-Champaign for RNA
sequencing. Samples were placed individually in a cryocentrifuge tube, labeled with a
roman numeral and sample number (Appendix I), packed in dry ice, and shipped via
FedEx.
52

e. RNA Sequencing:
RNA was sequenced using the Illumina Next Generation Sequencing platform,
sequencing by synthesis. This four-step process involves 1) preparation of an mRNA
library; 2) cluster amplification; 3) sequencing; and 4) alignment to a reference genome.
In step 1, a TruSeq Stranded mRNA library was prepared as outlined in the workflow
diagram in Figure 24. In overview, the total RNA for each sample was mRNA enriched
and then fragmented. A first strand of cDNA was synthesized followed by the synthesis
of the second strand of cDNA using dUTPs. The 3’ ends were adenylated and the
adaptors ligated. After ligation, the dUTPs were enzymatically removed to prepare for
PCR amplification. At the completion of the PCR amplification, the created library was
validated by the use of qPCR and Bioanalyzer. Lastly, the library was pooled and
normalized, resulting in the final library.

Figure 24: TruSeq Stranded mRNA library preparation workflow (Photo Credit: TruSeq
mRNA library Preparation Kit protocol). TruSeq stranded mRNA kit sequences RNA by
synthesis via Illumina Sequencing.

53

In step 2, the mRNA library was loaded into a flow cell and the fragments were
hybridized to the surface of the flow cell. The mRNA fragments were amplified into
clone clusters by bridge DE amplification. In step 3, sequencing by synthesis utilized
fluorescently-labeled nucleotides to sequence the RNA on the flow cell surface. For each
sequencing cycle, a single labeled dNTP was added to the nucleic acid chain thus serving
as a terminator for polymerization. Upon addition of each dNTP, the fluorescent dye is
imaged to allow for identification of the base and then is cleaved so the same process can
be performed on the next nucleotide. Since the base identity is determined one at a time
from signal intensity measurements during each cycle, the raw error rates are reduced,
resulting is highly accurate base-by-base sequencing (Illumina 2010). Finally, in step 4,
the data are aligned and compared to the reference genome as described below.
Upon completion of RNA sequencing, the number of reads per gene were
normalized due to potential differences in the number of reads as well as potential
differences in RNA composition. The Apis mellifera reference genome Amel_HAv.3.1
has a total of 12,090 genes. Due to low detectable expression, a total of 1,493 genes were
filtered out resulting in 10,597 genes that were eligible for differential expression
analysis (Valizadegan and Drnevich 2019).
f. Post-Sequencing Analysis and Statistics:
Results of RNA sequencing were reported by Negin Valizadegan and Jenny
Drnevich of the Keck Center for Comparative and Functional Genomics at University of
Illinois Urbana Champaign. The Apis mellifera transcriptome derived from genome
Amel_HAv3.1 and Annotation Release 104 from the National Center for Biotechnology
Information was used for quasi-mapping and gene counts (Valizadegan and Drnevich
2019). The Sequencing Center’s report was summarized using Multi-QC version 1.6
(Appendix J and K). A quality check on the raw data indicated that the reads were of high
quality and no adaptor sequences were found. Salmon version 0.13.1 was used to quasimap reads to the transcriptome and to determine the abundance of each transcript (Patro
et al. 2017). Approximately 70.1% to 82.3% of the reads were mapped to the
transcriptome. The unmapped reads were discarded (Valizadegan and Drnevich 2019). In
order to compare expression levels among the control and treatment groups, the number
of reads per gene were normalized using the EDE-R package, an essential step given that
54

there could be differences in the total number of reads and variations in RNA
composition.
To identify treatment effects, a multi-dimensional scaling plot was constructed
with the top 5,000 most variable genes using limma (Valizadegan and Drnevich 2019).
Differential gene expression analysis was conducted using the limma-trend method
(Valizadegan and Drnevich 2019) by comparing the following pairs: Treatment-0h versus
Treatment-4h, Control-0h versus Control-4h, Treatment-0h versus Control-0h, and
Treatment-4h versus Control-4h. A test for any interaction between treatment and time
was also performed. To adjust for multiple testing, a “global” False Discovery Rate
(FDR) correction was done across p-values for all comparisons. The number of upregulated and down-regulated genes for each pairwise comparison was made using an
FDR p-value of <0.1. A heatmap was made to visualize the overall gene expression
patterns with a significant one-way ANOVA, FDR <0.05.
A total of 10,597 genes were returned post-sequencing. Of the 10,597 genes, a
total of 7,806 genes were considered characterized (ie., functional and/or pathway
information available and/or has a mammalian homolog). The gene expression
differences among the control and treatment groups were determined using R
programming with the use of one-way ANOVA corrected with Bayesian probabilities. In
our laboratory, an automated algorithm was developed to identify genes where the
expression data was inconsistent for more than two out of the six bee samples within a
specific pairwise comparison. Any genes where >33% of the data was inconsistent (over/under-expression among the samples) were removed from that treatment group. Table 2
breaks down the treatment comparison and the number of genes removed from each
comparison. The specific outlier genes removed from each comparison can be found in
Appendix L. The automated algorithm was written in the Python 2.7 programming
language (Van Rossum and Drake 1995, Appendix M).

55

Table 2: Number of genes removed for each comparison. A total of 10,597 genes were identified after
RNA sequencing. A bee sample was deemed an outlier if it differed greatly from the other bees in the
assigned treatment group. A gene was removed from analysis if more than two bees out of the six bee
samples for a specific gene were identified as going in the opposite direction.

Treatment Group

Number of Genes
Removed
15

Treatment-4h versus
Treatment-0h
Control-4h versus
Control-0h
Treatment-0h versus
Control-0h
Treatment-4h versus
Control-4h

80
25
137

g. Gene Enrichment and Pathway Analysis:
Changes in gene expression among control bees and bees exposed to 0.9 ng/bee
(1/20th of LD50) of imidacloprid were examined by making four pairwise comparisons
(Figure 25).

Figure 25: Pairwise comparisons of control and imidacloprid-treated bees. The schematic illustrates
the statistical comparisons made among four treatment groups.

56

The groups Control-0h and Control-4h were compared to assess research design
effects while Control-0h and Treatment-0h were compared to asses randomization and
the small sample size. Therefore, if a gene is shown to be up-regulated or down-regulated
in one of these comparisons, it is not due to the treatment and is an artifact of the
experimental design. Our main focus is on the effects of pesticide treatment, that is, in the
differential gene expression of Treatment-4h compared to Treatment-0h and Control-4h
compared to Treatment-4h.
DAVID, standing for Database for Annotation, Visualization, and Integrated
Discovery, was utilized for Functional Analysis (Nature Protocols 2009). DAVID
Bioinformatics Resource 6.8 uses gene set enrichment to provide functional annotation,
gene functional classification, gene ID conversion, and gene name batch viewer for
individuals to use in order to understand biological meaning behind a large list of genes.
The identification numbers for were entered into the Gene List, identified as Entrez ID
numbers, determined to be a gene list, and submitted. This process was performed for all
significant (p<0.1) down-regulated and up-regulated genes separately. A function or
functional cluster that had an FDR p<0.05 was identified as significant. The function of
the cluster, group enrichment score, and FDR were documented. The overall group
enrichment score is based on the modified Fisher’s Exact p-value for each term member.
The higher the enrichment score, the more enriched, and the smaller the p-value (Nature
Protocols 2009). This process was performed for Treatment-4h versus Treatment-0h and
Control-4h versus Control-0h.
Standing for Kyoto Encyclopedia of Genes and Genomes, KEGG Mapper (KEGG
Mapper 2020) was used for Pathway Analysis. This database is a resource that helps to
visualize the interactions among genes in various pathways and shows how a change in
the function of gene can alter the function of genes further downstream. To retrieve a list
of pathways, a gene list was inputted into KEGG. For this research, all 10,597 genes were
entered. If a gene was considered significant (p<0.1), the gene was designated purple
(down-regulated) or pink (up-regulated). If the gene was not significant (p>0.1) then it
was designated orange (up-regulated) or red (down-regulated). After the genes were
entered, the analysis searched the knowledge base and identified pathways that contained
a gene from the list that was inputted. The genes, both significant and not significant,
57

appeared in the pathway with the appropriate color designation. A pathway was deemed
interesting based on previous studies and literature about the effects of neonicotinoids on
the honey bee. This process was performed for Treatment-4h versus Treatment-0h and
Control-4h versus Control-0h.

IV.

Results

Experiment 1: The Effects of Sublethal Doses of the Neonicotinoid Pesticide Imidacloprid
on Motor Function and Cellular Responses in Honeybees
Motor Function Responses of Honeybees:
Figure 26 shows the motor responses of honeybees exposed to sublethal doses of
imidacloprid ranging from 1/5th to 1/500th of the LD50. A motor score of 8, represented
normal leg and abdomen movement, antennae responsiveness, and a proboscis extension
reflex. No differences were observed in the motor scores of the bees in the positive (heat
shock) and negative (room temperature) control groups (p<0.05). Motor scores in both
positive and negative controls were normal, near 8. Lower motor scores, less than 8, were
observed in honeybees that were treated with sublethal doses of imidacloprid. All of the
scored motor functions were impaired relatively to the same degree. Honeybees in the
treatment groups 1.8 ng/bee (1/10th LD50) and 3.6 ng/bee (1/5th LD50) had lower motor
scores than negative and positive controls (p <0.05). Motor responses of the bees exposed
to imidacloprid doses ranging from 0.9 ng/bee (1/20th of the LD50) to 0.036 ng/bee
(1/500th of the LD50) did not differ compared to controls.

58

Figure 26: Average motor scores of honeybees exposed to sublethal doses of imidacloprid. Honeybees
were exposed to serial dilutions of the oral LD50 of imidacloprid (18.0 ng/bee) or 1.5 M sucrose (controls).
Bars representing means ± standard error (n=261). Chi-square analysis was performed (SPSS Statistics).
Means connected by the blue line did not differ from the negative control. Means connected by red line did
not differ from the positive control. Means connected to the black line did not differ from the 1/100th LD50
treatment (SOD only). Broken lines with asterisk(s) indicate where means differed from a comparison:
*=P<0.05, **=P<0.01. NS represents no significant differences.

Dose-response analysis of the percentage of bees that were impaired (motor
scores < 8) per treatment showed a strong sigmoidal curve fit with an estimated toxic
dose at which 50% of bees were impaired (TD50) of 1.78-2.25 ng/bee (Figure 27). A
supplemental table of the average motor scores is provided in Appendix N.

59

Figure 27: Imidacloprid motor response curve. Percentage of bees that were impaired were plotted at
each sublethal dose of imidacloprid to determine the toxic dose 50 (TD50). A nonlinear regression fitted a
sigmoidal dose-response curve in GraphPad v6 by Prism.

Superoxide Dismutase (SOD) Activity in Honeybees:
Figure 28 shows the SOD activity in homogenized honeybee head capsules
exposed to sublethal doses of imidacloprid ranging from 1/5th to 1/500th of the LD50. No
differences were observed in the SOD activity in the bees in the positive (heat shock) and
negative controls (room temperature). The results of the Tukey-Kramer post-hoc test
indicated elevated SOD activity in both the 1/5th (3.6 ng/bee) and the 1/10th (1.8 ng/bee)
doses when compared to the 1/100th (0.18 ng/bee) dose (ANOVA: F=3.2897, p=0.0037).

60

Figure 28: SOD activity in honeybees treated with imidacloprid. There is a significant difference in
the SOD activity between the 1/5th dose (3.6 ng/bee) and the 1/100th dose (0.18 ng/bee) and between the
1/10th dose and the 1/100th dose of imidacloprid (ANOVA: F=3.2897, p=0.0037). Means connected by
the blue line did not differ from the negative control. Means connected by red line did not differ from
the positive control. Means connected to the black line did not differ from the 1/100th LD50 treatment
(SOD only). Broken lines with asterisk(s) indicate where means differed from a comparison: *=P<0.05.

HSP70 Concentrations in Honeybees:
The levels of the cellular stress marker HSP70 in bees exposed to sublethal doses
of imidacloprid are shown in Figure 29. The positive controls, bees that were heat
shocked at 42°C, had increased levels of HSP70 compared to the negative controls
(samples that were fed 1.5 M sucrose and left at room temperature). The Tukey-Kramer
post-hoc test indicated that HSP70 levels in the 0.18 ng/bee (1/100th dose of LD50) and
3.6 ng/bee (1/5th dose of LD50) treatments were lower (p<0.05) than the positive control.
The 3.6 ng/bee (1/5th dose of LD50) treatment disrupts the cellular stress response of the
honeybees. Doses of 1.8 ng/bee (1/10th of LD50) to 0.36 ng/bee (1/50th of LD50) showed
intermediate levels of stress that overlapped both positive and negative controls.

61

Figure 29: The effects of sublethal doses of imidacloprid on HSP70 levels in honeybees. One-way
ANOVA (JMP Pro 12.2.0) indicated significant differences among groups. Means connected by the blue
line did not differ from the negative control. Means connected by red line did not differ from the positive
control. Means connected to the black line did not differ from the 1/100th LD50 treatment (SOD
only). Broken lines with asterisk(s) indicate where means differed from a comparison: ***=P<0.001.

Overall Results for Experiment 1:
Data from all of these tests, motor response scores (Figure 26 and 27), SOD
activity (Figure 28), and HSP70 concentrations (Figure 29) indicate that a dose of 0.9
ng/bee was intermediate among doses with an elevated stress response that was not
associated with impaired motor coordination. Therefore, a conservative sublethal
imidacloprid dose of 0.9 ng/bee (1/20th of LD50) was selected to treat bees for RNA
transcriptome analysis for Experiment 2. We hypothesized that an intermediate dose of
0.9 ng imidacloprid/bee would induce a sublethal stress response but not likely elicit a
massive apoptotic response, typical of an imidacloprid exposure from which bees might
recover from in nature.
Experiment 2: The Effect of Sublethal Doses of Imidacloprid on Gene Expression in the
Brain
Gene Expression Analysis: RNA-Sequencing:
The RNA sequencing analysis report, shown in Appendix J, contains the output
from the raw read files. FastQC (Appendix K) was used to generate quality check reports
for the sequencing reads of all 24 honeybees (n=6 per treatment group). The average per62

base read quality scores for all samples had a Phred score, a measure of the quality of the
identification of the nucleobases generated by automated DNA sequencing, over 30. This
value indicates that all of the samples after sequencing were of high quality, with a less
than 1 in 1,000 probability of an incorrect base call and a 99.9% base call accuracy
(Technical Note 2011; Appendix K).
Gene expression patterns in the four groups were analyzed by determining
differential gene expression (DGE) in pairwise comparisons. We predicted that the group
treated with imidacloprid, the Treatment-4h group, would have the greatest number of
differentially expressed genes compared to the other three groups, Control-0h, Control4h, and Treatment-0h. The DGE results however, do not match our prediction. Very few
genes were differentially expressed in the Treatment-4h group compared to the
Treatment-0h group (121 genes, sum of down-regulated and up-regulated genes).
Surprisingly, there were many more differentially expressed genes in pairwise
comparisons involving the Control-4h group than comparisons between the other
treatment groups. For example, the number of differentially expressed genes were greater
in Control-4h group compared to the Control-0h (4,285 genes, sum of down-regulated
and up-regulated genes) groups and the Treatment-4h and Control-4h (3,872 genes, sum
of down-regulated and up-regulated genes). Few genes were differentially expressed in a
comparison of Treatment-0h and Control-0h (232 genes, sum of down-regulated and upregulated genes). This result was anticipated because these groups, the Control-0h and
Treatment-0h, were included in the experimental design to assess randomization and the
small sample size. Therefore, the 232 genes that are up-regulated or down-regulated in
this comparison, are not due to the treatment but rather are artifacts of the experimental
design.
Differential gene expression (DGE) was performed for all four pairwise
comparisons including a test for interaction between treatment and time (Table 3)
(Valizadegan and Drnevich 2019).

63

Table 3: Differential Gene Expression (DE) for all four pairwise comparisons (Valizadegan and
Drnevich 2019). The number of down-regulated, up-regulated, and non-significant genes were identified
among all four pairwise comparison. The total number of significant differentially expressed genes is the
sum of the significant down-regulated and up-regulated genes in each comparison (p<0.1).

Down-Regulated
Not Significant
Up-Regulated
Total Number of
Significant
Differentially
Expressed Genes

T4vsT0
46
10,476
75
121

C4vsC0
2,369
6,312
1,916
4,285

T0vsC0
35
10,365
197
232

T4vsC4
1,663
6,725
2,209
3,872

Interaction
1,091
7,805
1,701
2,792

There are many more differentially expressed genes in pairwise comparisons involving
the Control-4h group than comparisons between the other treatment groups. Like the
DGE data, the multi-dimensional scaling plot clearly shows that gene expression patterns
in the Control-4h group differs from the other three groups. The Control-4h group is
sharply separated on the X axis from the other three groups. To illustrate overall
expression patterns, a one-way ANOVA was calculated across all four groups to select
genes to visualize in a heat map. Using an FDR <0.05, a total of 3,819 genes were
identified as different across the four treatment groups (Figure 31). The comparisons that
are most important to examine are Control-4h versus Control-0h and Treatment-4h versus
Control-4h. Control-4h versus Control-0h comparison reveals the effects of experimental
manipulation on control animals; therefore, there should not be a large change in
differential expression in comparison to the treatment group (Treatment-4h versus
Control-4h). The heat map reveals that the Control-4h shows a greater change in gene
expression when compared to Control-0h, Treatment-0h, and Treatment-4h.
To adjust for multiple testing that occurred due to the increased number of
differentially expressed genes in the Control-4h, a “global” False Discovery Rate (FDR)
correction was performed across the p-values for all five comparisons. This ensured that a
gene with the same raw p-value in two different comparisons would not end up with
extremely different FDR p-values.

64

Figure 30: One-way ANOVA heat map of Control-0h, Control-4h, Treatment-0h, and Treatment-4h
(Valizadegan and Drnevich 2019). Each lane within a group represents an individual bee (1-6). Genes
that have been down-regulated are colored blue, up-regulated genes are colored red, and white denotes no
change in gene expression. The color intensity represents the degree of gene expression change.

65

Multi-dimensional scaling, a method used to visualize the level of similarity of
individual cases of a dataset at a higher level, was used to identify possible treatment
effects. The multi-dimensional scaling plot in Figure 31 is based on 5,000 of the most
variable genes. Dimension 1 (X-axis) illustrates approximately 25% of the total
variability and explains the separation of the Control at Time 4h samples from the other
groups. Dimension 2 (Y-axis) encompasses approximately 17% of the variability. This
does not explain groupings of the effect of treatment, control, and/or the time, however
(Valizadegan and Drnevich 2019). It would be expected for the multi-dimensional scaling
plot to show Control-0h, Control-4h, and Treatment-0h to be clustered relatively close
while Treatment-4h showed greater separation. The overall findings of this plot;
however, show that Control-4h specimens (depicted by the red circle) had the greatest
variability in comparison to Control-0h, Treatment-0h, and Treatment-4h.

Figure 31: Multi-dimensional scaling plot performed on 5,000 of the most variable genes
(Valizadegan and Drnevich 2019). Each individual bee in the four groups is represented on the plot (n=6
for each group). Each of the four groups is represented by a different color. Control-0h bees are labeled
Neg_0 followed by the sample bee number and appear in black. Control-4h bees are labeled Neg_4
followed by the sample bee number and appear in red. Treatment-0h hour bees are indicated by Treat_0
followed by the sample bee number and are in green. Treatment-4h bees are indicated by Treat_4 followed
by the sample bee number in blue.

66

A list of the top 100 most significant differentially expressed genes for Treatment4h and Control-0h can be found in Appendix O. All 10,597 honeybee genes were
analyzed in KEGG. The Control-4h versus Control-0h is of interest since this comparison
contained the greatest differential gene expression. The schematic from Figure 25
illustrates that this should show effects from research design; however, further
investigation into the large differential gene expression difference is needed. Table 4 lists
key genes and pathways that were significantly changed in the Control-4h versus
Control-0h group comparison.
Table 4: Key Pathways and genes at Control-4h and Control-0h. Pathway analysis was performed
using KEGG. Significant genes had a p-value <0.1.

Key KEGG
Pathways
Circadian
Rhythm-Fly
Metabolism of
Xenobiotics by
Cytochrome
P450
Drug
MetabolismCytochrome
P450
Glutathione
Metabolism
FOXO
Signaling
Pathway

Peroxisome

Longevity
Regulating
PathwayMultiple
Pathways

Apoptosis-Fly

Upregulated

KEGG Genes
Period circadian protein (Per)

DownRegulated

KO #

√

K02633

Glutathione s-transferase D1 (GstD1)
Glutathione s-transferase (GstS1)

√
√

K00799
K04097

Glutathione s-transferase D1 (GstD1)
Glutathione s-transferase (GstS1)

√
√

K00799
K04097

Glutathione s-transferase D1 (GstD1)
Glutathione s-transferase (GstS1)
Catalase (Cat)
Insulin-like receptor-like (InR-2)
Phosphatase and tensin-like (Pten)
Superoxide dismutase 2, mitochondrial
(Sod2)
Catalase (Cat)
Dehydrogenase/reductase SDR family
member 4 (DHRS4)
Fatty acyl-CoA reductase 1 (FAR1)
Superoxide dismutase 1 (Sod1)
Superoxide dismutase 2, mitochondrial
(Sod2)
Adenylate cyclase 3 (Ac3)
Catalase (Cat)
Heat shock protein cognate 4 (Hsc70-4)
Insulin-like receptor-like (InR-2)
Superoxide dismutase 1 (Sod1)
Superoxide dismutase 2, mitochondrial
(Sod2)
Broad-complex (Br-c)
Cytochrome c (CytC)
Ecdysteroid-regulated gene E74 (E74)
Ecdysone receptor isoform A (Ecr)

√
√

K00799
K04097
K03781
K04527
K01110
K04564

67

√

√
√
√

√
√
√
√

√
√
√
√
√
√

√
√
√

√
√

K03781
K11147
K13356
K04565
K04564

K08043
K03781
K03283
K04527
K04565
K04564
K02174
K08738
K09428
K14034

Oxidative
Phosphorylation

Mushroom body large-type Kenyon cellspecific protein 1 (Mblk-1)
Ultraspiracle (Usp)
ATP synthase lipid-binding protein,
mitochondrial (ATP5G2)
Cytochrome c oxidase subunit VIb
polypeptide 1 (Cox6b1)
Cytochrome c oxidase subunit 6c (Cox6bc)
NADH dehydrogenase (ubiquinone) 1 beta
subcomplex, 2, 8kDa (Ndufb2)
NADH-ubiquinone oxidoreductase 75 Dka
subunit, mitochondrial (Ndufs1)
NADH dehydrogenase [ubiquinone] ironsulfur protein 5 (Ndufs5)
Cytochrome b-c1 complex subunit 10
(Uqcr11)
Vacuolar H+ ATP synthase 16 Dka
proteolipid subunit (Vha16)

√

K20015

√
√
√
√
√

K14030
K02128
K02267
K02268
K03958
√

√
√

K03934
K03938

√

K00420
K02155

Significant genes with a p<0.1 in the comparison Control-4h versus Control-0h
and Treatment-4h versus Treatment-0h were analyzed with DAVID. Just as the KEGG
analysis showed, the Control-4h versus Control-0h comparison contained the greatest
amount of gene enrichment. Although both DAVID and KEGG analysis yielded similar
findings, these findings are inconsistent with the expectation that Treatment-4h versus
Treatment-0h would reflect the greatest change due to pesticide treatment. The results of
the DAVID analysis showed that Zinc Finger, Transmembrane, GTP-binding, Ubiquitin
Protein Transferase Activity, ATP-binding, and Intracellular Signal Transduction had
significant down-regulated enrichment changes in genes in the Control-4h versus
Control-0h comparison. DAVID analysis also found that Ribosome, Small Nucleolar
Ribonucleoprotein, Large Ribosomal Translation, Nucleotide Binding, Immunoglobulin,
RNA Polymerase, DNA-templated Transcription, and Pheromone/General Odorant
Binding had significant up-regulated enrichment changes in genes in the Control-4h
versus Control-0h comparison.

V.

Discussion

Motor Function Responses:
Bees exposed to sublethal doses of imidacloprid showed impaired motor
responses compared to positive and negative control bees at doses of 1.8 ng/bee (1/10th
LD50) and 3.6 ng/bee (1/5th LD50). Motor responses did not differ between the positive
and negative control groups. This was expected since the positive control group, while
68

subjected to heat stress, was incubated at a temperature about 1°C below the critical
thermal maximum for honeybees (Atmowidjojo et al. 1997). Motor responses were still
intact at this temperature. Sigmoidal dose-response model predicted that 50% of the
forager population had impaired motor responses between 1.78-2.25 ng/bee of
imidacloprid (Figure 27). This range of values includes the 1/10th LD50 of imidacloprid
(1.80 ng/bee) in our experiment.
Because of their widespread use, water solubility, persistence in the soil, and
uptake by plants, neonicotinoid pesticides are broadly present in the environment.
Neonicotinoids have been found in trace levels in plants, pollen and nectar and have been
detected in bees wax, honey, and honeybees themselves (reviewed in Blacquière et al.
2012). While the levels of neonicotinoids employed in agriculture are not lethal to bees,
chronic sublethal exposures may occur. For example, honeybees can reach our observed
TD50 value of 1.78 to 2.25 ng in a few trips to common field plants treated with
imidacloprid. Stoner et al. (2012) found that squash nectar contains 10 μg/L, or 10 pg/μL
of imidacloprid. Sylvester et al. (1983) reported that bees can carry up to 62 μL of nectar.
This would mean bees can carry 0.620 ng of imidacloprid from a full trip to a squash
plant. Based on this, in less than three trips worker honeybees could experience the TD50
value of 1.78 to 2.25 ng/bee. Girolami et al. (2009) found that corn guttation holds up to
200 mg/L, or 200 ng/μL of imidacloprid. One exposure to corn guttation would be lethal
and partial consumption of only 0.1 μL would lead to dosages greater than the TD50 value
of 1.78-2.25 ng/bee. These sublethal exposures combined with other stressors fit the
criteria of the Sublethal Stress Model (Bryden 2013) and can impair the health of the
hive.
The declines in motor responses of honeybees exposed to sublethal doses of
imidacloprid may have significant consequences. Flight patterns, behavior, and foraging
ability have all been shown to be reduced in honeybees exposed to imidacloprid (Tan et
al. 2014, Karahan et al. 2015, Decourtye et al. 2004a, Decourtye et al. 2004b). This
could be attributed to the impaired motor function of honeybees. Worker bees who do not
forage efficiently will take longer to return to the hive and worker bees that cannot fly
properly may not return to the hive at all. Impairments in the proboscis extension reflex
and in antennae movement could lead to problems locating and foraging nectar. The
69

reduction in antennae function could diminish the ability of the bees to detect location,
resulting in altered flight patterns.
The overall impairment of worker bees can negatively affect the health of the hive
by lowering its nectar and pollen supplies. Contaminated food in the hive can also lead to
weakened larvae that will grow into the next generation of worker bees. This may also
lower the worker bees’ ability to defend the hive. Guard bees specialize in olfactionbased nestmate recognition and alarm-pheromone-mediated recruitment of nestmates to
sting. Stinging is determined by visual, tactile, and olfactory stimuli (Hunt 2007).
Overall, the impaired motor functions of worker bees are detrimental to a honeybee
colony. A decline in motor responses due to imidacloprid intoxication could result in
delayed decision-making and other cognitive skills (Tan et al. 2014). When these factors
are compounded with other sublethal stressors such as mites, nutritional stress, or
environmental stress, it is evident that the effects on the hive may be devastating.
Superoxide Dismutase (SOD) Activity:
SOD is a potent antioxidant enzyme that is a marker of oxidative stress. SOD
combats oxidative stress by scavenging the superoxide anion, a reactive oxygen species,
and converting it to hydrogen peroxide, a substrate of catalase. Elevated levels of SOD
have been proposed to serve as an organism’s first line of defense to prevent the
damaging effects of superoxide on cellular biomolecules (Ighodaro and Akinloye 2018).
In our study, SOD activity was higher in bees treated with the 1/5th (3.6 ng/bee) and the
1/10th (1.80 ng/bee) doses of imidacloprid compared to those that received the 1/100th
(0.180 ng/bee) dose (Figure 28). As seen in the case of motor responses, no significant
difference was observed in SOD activity between the positive (incubated at 4 hours at
42°C) and negative control (room temperature). This response is to be expected since
SOD is not affected by temperature.
Several studies reported an increase in the activities of oxidative enzymes such as
catalase, glutathione-S-transferases, superoxide dismutase, and xanthine oxidase in bees
after sublethal exposure to pesticides such as pyrethroids and organophosphates
(Chakrabarti et al. 2015). Chakrabarti et al. (2015) observed in field studies of native
bees that SOD levels were greater in bees from agricultural areas with high pesticide use
compared to those from pesticide-free zones. Oxidative stress is known to influence the
70

survival and homing ability of honeybees; increases in oxidative stress resulted in
decreased survival and foraging abilities of the affected honeybees (Simone-Finstrom et
al. 2016). Our findings, elevated SOD levels following exposure to sublethal doses of
imidacloprid, are indicative of oxidative stress and suggest that honeybee foragers and
the overall health of the colony could be negatively impacted.
Heat Shock Proteins (HSP70):
Heat shock proteins function as molecular chaperones to prevent misfolding of
cellular proteins during conditions of heat stress or other types of cellular stress,
including exposure to toxins (Feder et al. 1997). The positive controls, bees that were
heat shocked at 42°C for four hours, had increased levels of HSP70 compared to the
negative controls (samples that were fed 1.5 M sucrose and remained at room
temperature). HSP70 levels in the 0.18 ng/bee (1/100th dose of LD50) and 3.6 ng/bee
(1/5th dose of LD50) treatments were lower than the positive control (ANOVA values, p
<0.05). The 3.6 ng/bee (1/5th dose of LD50) treatment had a decreased cellular stress
response as reflected by HSP70 levels while doses of 1.8 ng/bee (1/10th of LD50) and 0.3
ng/bee (1/50th of LD50) showed intermediate levels of stress that overlapped both positive
and negative controls. These findings support the hormesis model with a doserelationship that is stimulatory at low doses and inhibitory at higher doses (see Figures 11
and 29). The 1/500th dose appears to be an exception. We cannot account for the
discrepancy but this could be due to errors when making up the solution or other
unknown regulatory mechanisms.
Previous studies have shown that HSP70 is a useful predictor of cellular stress in
honeybees. Chacon-Almeida et al. (2000) observed an induction of heat shock proteins in
the larval fat body of honeybees incubated at 42°C for 2 hours. Work performed by
Hranitz and colleagues (2010) demonstrated that ethanol doses that affect honeybee
learning and behavior cause significant increases in HSP70 concentrations, a cellular
stress response. The dose-response relationship between ethanol concentration and
HSP70 levels also followed the pattern of the hormesis model. In a study by Sahebzadeh
and Lau (2017), worker bees infested with Varroa mites or exposed to sublethal
acaracides such as thymol or eucolyptol showed a dose dependent up-regulation of all
HSPs, pointing to their utility as biomarkers when bees are exposed to toxic or
71

pathogenic stress. These studies lend support to our findings that the elevated levels of
HSP70 that we observed in honeybees exposed to sublethal doses of imidacloprid are
indicative of cellular stress.
Determination of Imidacloprid Dose to Use for Experiment 2:
Sublethal doses of imidacloprid impaired motor responses and produced elevated
levels of HSP70 (marker of cellular stress) and SOD (marker of oxidative stress) in
honeybees (Figures 26-29). As mentioned in the Introduction, the hormesis model
illustrates parabolic curves with positive, negative, and peak response levels across a
spectrum of doses. Lower doses of imidacloprid elicited a beneficial, protective effect
within the hormetic zone while a 3.6 ng/bee (1/5th dose of LD50), a high dose, resulted in
inhibition. A conservative sublethal imidacloprid dose of 0.9 ng/bee (1/20th of LD50) was
selected to treat bees for RNA transcriptome analysis for Experiment 2. This dose of 0.9
ng/bee was intermediate among doses with an elevated stress response. We estimated that
0.9 ng imidacloprid/bee would induce a sublethal stress response but not likely produce a
massive apoptotic response, typical of an imidacloprid exposure from which bees might
recover from in nature.
Mortality:
The experiment-wide mortality rate of the restrained bees after 22-24 hours for
Experiment 1 was 38%. Despite the mortality rates, negative and positive controls were
significantly different in stress response. There was 0% mortality of bees in Experiment 2
treatments.
RNA-Seq Results:
Collectively, the DE analysis, Multi-dimensional Scaling Plot, and the one-way
ANOVA heat map agree that the Control-4h group had the greatest changes in gene
expression compared to the other groups. This is counter to our prediction that the
greatest changes would be in the bees treated with imidacloprid. Our Control-4h samples
showed elevations in SOD gene expression (Table 3). This finding contradicts our
independent results from Experiment 1. Specifically, Experiment 1 showed that SOD
activity was elevated in bees treated with imidacloprid compared to controls.

72

The increased gene expression in the Control-4h group that was observed is also
inconsistent with previous studies. For example, a study of similar design by the Hranitz
laboratory showed that RNA expression was affected by ethanol. A total of 609
microarray ESTs in the honeybee brain displayed different expression to acute (4 hour)
ethanol treatment (personal communication). Furthermore, a review of the imidacloprid
literature shows varied and significant impacts of this pesticide on the physiology and
behavior of honeybees that would likely be reflected in genome-wide responses (Henry et
al. 2012, DiPrisco et al. 2013, Wu-Smart et al. 2012, Pereira et al. 2020, and Decourtye
et al. 2004). Recent studies investigating the effects of neonicotinoids on the honeybee
transcriptome have also shown significant changes in gene expression in bees following
neonicotinoid treatment. For example, Shi et al. (2017) reported a total of 609
differentially expressed genes when bees were treated with the second-generation
neonicotinoid thiamethoxam.
Post-Hoc Hypothesis and Supporting Evidence:
There is a wealth of corroborating evidence that suggests that the Control-4h and
the Treatment-4h samples were mislabeled and switched. These discrepancies have led us
to propose a post-hoc hypothesis that the Control-4h and Treatment-4h samples were
mislabeled prior to shipping to University of Illinois Urbana-Champaign. RNA isolation
of Control-4h and Treatment-4h bees were performed on the same day so this mishap
may have occurred at that time.
Several steps were taken to investigate the possibility of a switch between the
Control-4h and Treatment-4h groups. First, the laboratory notebook was checked to
confirm that the sample numbers and RNA concentration values matched the document
that was sent with the specimens to University of Illinois Urbana-Champaign. Second,
the Control-4h and Treatment-4h specimens that were not sent for RNA sequencing were
re-analyzed for RNA content. The repeated RNA yields were not consistent with the
original measurements; however, the RNA concentrations were consistent between the
two measurements of the same samples about one month apart. Third, University of
Illinois Urbana-Champaign repeated RNA sequencing on the remainder of each original
bee sample. As a final attempt to understand the results from RNA sequencing, further
data analysis was performed on the output provided by University of Illinois Urbana73

Champaign. A Python code was used to identify outliers within the same treatment group
for a given expression level. A gene or honeybee was deemed an outlier if it met one of
the following criteria: if more than two bee samples out of the six samples for a specific
gene were identified to be going in the opposite direction (positive and negative) or if the
read value returned was greater than 2 standard deviations (SD) different from the other
honeybee samples with that treatment group. The results of this data cleaning with
outliers removed can be found in Table 5.
Table 5: Removal of genes and outliers using Python. The number of genes removed from each
treatment comparison and the number of outliers or honeybees removed per treatment comparison from the
gene list is reflected below. Significance was p<0.1.

Number of
Genes
Removed
Number of
Outliers

T4vsT0
15

T0vsC0
25

C4vsC0
80

T4vsC4
137

2,669

2,804

2,804

2,669

The results in Table 6 show that the removal of outliers from the dataset did not change
the interpretation of the differential gene expression. The Control-4h group still displayed
most of the changes in gene expression. Based on the supporting evidence, we feel
confident that the Control-4h and Treatment-4h samples were switched and mislabeled.
The pair-wise comparisons of the Control-4h and Treatment-4h bees have been renamed
to reflect this mislabel (Table 6). This label change will remain for the rest of the
discussion. The gene expression data in Table 6 has been cleaned by the removal of
outliers in Table 5.

74

Table 6: Results of further data cleaning using Python. The total number of genes that are up-regulated
and down-regulated are based on the genes that are considered significant. The number of genes for each
comparison takes into account the outliers that were removed from the gene list. Significance was based on
p-value <0.1.

DownRegulated
Not
Significant
UpRegulated
Genes
Removed

C4vsT0
37

T0vsC0
32

T4vsC0
2,353

T4vsC4
1,569

10,476

10,365

6,312

6,725

69

175

1,852

2,166

15

25

80

137

Pathway and Gene Set Enrichment Functional Analysis:
The molecular responses by the honeybee brain transcriptome to acute sublethal
doses of imidacloprid can be interpreted by examining statistical comparisons among the
four groups. Our focus is on comparing the differential expression between Treatment-4h
and Control-0h, post label change, which directly looks at the effects of pesticide
treatment. Of the 10,597 total genes, 4,205 genes were considered significant (2,353
down-regulated genes and 1,852 up-regulated genes). Key cellular pathways that were
affected by imidacloprid treatment were the peroxisome pathway, metabolism of
xenobiotics by cytochrome P450, circadian rhythm, drug metabolism, FOXO signaling,
longevity regulating pathways, apoptosis, and oxidative phosphorylation (Table 4).
The peroxisome pathway provides a clear example of the cellular stress produced
by the exposure of bees to sublethal doses of imidacloprid (Figure 32). In the Control-4hTreatment-0h comparison (Appendix P), both catalase and SOD were down-regulated,
indicative of little to no cellular stress. However, in the Treatment-4h-Control-0h
comparison, catalase is down-regulated while SOD is up-regulated. This is confirmed by
the SOD assay performed in Experiment 1. This provides additional supporting evidence
that a switch occurred between Treatment- 4h and Control-4h. If a switch did not occur,
then SOD would have been elevated in the controls which is opposite of what was shown
in the SOD assay.

75

Figure 32: Peroxisome Pathway of Treatment-4h-Control-0h. KEGG was used as the enrichment
analysis database. Genes colored red are considered down-regulated (not significant) and genes colored
orange are considered up-regulated (not significant) (significance p<0.1). Enzymes are indicated by the
green color.

Seven genes involved in circadian rhythm are found in the honeybee genome.
One of the major ones, the Period (per) gene was dramatically affected when exposed to
76

imidacloprid. A comparison between Treatment-4h-Control-0h revealed significant
down-regulation of the Per gene (Figure 33). The Control-4h-Treatment-0h comparison
(Appendix P) showed a non-significant down-regulation of Per. Changes in circadian
locomotor rhythms, mating behavior, and development time from egg to adult are
affected by mutations in Per (Toma et al. 2000). While mRNA levels of Per were found
to shift in bees of all ages, the most significant changes could be found in foragers (Toma
et al. 2000). Karahen et al. (2015) found that imidacloprid reduced foraging activity by
honeybees, where bees spent more time between trips in the hive. Tasman et al. (2020)
investigated the effects of imidacloprid on bumble bees and discovered that
concentrations as low as 1μg/L caused reductions in locomotion and foraging activities
and mistiming that resulted in increased foraging at night and increased daytime sleeping.
These changes in honeybee behavior could have dramatic effects on the survival of the
colony.

Figure 33: Circadian Rhythm Pathway-Fly Treatment-4h-Control-0h. KEGG was used as the
enrichment analysis database. Genes that are colored purple are considered significantly down-regulated
(significance p<0.1). Enzymes are indicated by the green color.

Forkhead box transcription factors (FOXO) are proteins involved in a wide range
of cellular functions such as cellular differentiation, apoptosis, cell cycle regulation,
autophagy, DNA damage and repair, and as mediators of oxidative stress (Farhan et al.
77

2017). The Treatment-4h and Control-0h comparison revealed a significant reduction in
expression of the PTEN gene (Figure 34). An inhibition of PTEN enzyme activity (acts
as a tumor suppressor by regulating cell division) leads to a decrease in FOXO. During
oxidative stress, FOXO induces autophagy, a process which clears accumulated toxins
and promotes cell survival (Maiese 2015, U.S. Library of Medicine 2015). In the case of
imidacloprid-treated bees where FOXO is reduced, there may be less autophagy,
resulting in a greater accumulation of toxins and lowered cell survival. Oxidative stress
has also been reported to modify the interactions of FOXO with other proteins that
ultimately can affect the survival of neurons (Maiese 2015). The interpretation of the
FOXO pathways is further complicated by its acetylation and numerous different forms.

Figure 34: FOXO Signaling Pathway Treatment-4h-Control-0h. KEGG was used as the enrichment
analysis database. Genes that are colored red are considered down-regulated (not significant), genes that are
colored purple are significantly down-regulated, and genes that are colored orange are considered upregulated (not significant) (significance p<0.1). Enzymes are indicated by the green color.

Detoxification enzymes are essential for honeybees to remove xenobiotics, such
as neonicotinoids, that they are exposed to in their environment. Possessing roughly 3050% fewer detoxification genes than Anopheles and Drosophila, bees are particularly
vulnerable to the effects of xenobiotics. Xenobiotics are detoxified by three superfamilies
of enzymes carboxylesterase (CCE), cytochrome P450 (P450), and glutathione Stransferase (GST) (Honeybee Genome Sequencing Consortium 2006, du Rand et al.
78

2015). As seen in Table 4, pathways for drug metabolism-cytochrome P450, metabolism
of xenobiotics by cytochrome P450, and glutathione metabolism are significantly upregulated in the imidacloprid treated group (0.9 ng/bee of imidacloprid) compared to the
control group. This was due to the up-regulation of glutathione-S-transferases, a family
of Phase II detoxification enzymes that catalyze the conjugation of glutathione to a
variety of toxic compounds for the purpose of detoxification. In addition to its role in
detoxification, glutathione-S-transferases have been shown to protect bees against
oxidative damage caused by ROS (Yan et al. 2013). Glutathione-S-transferase has
previously been shown to be up-regulated in honeybees following exposure to nicotine
(du Rand et al. 2015). The up-regulation of the xenobiotic, drug metabolism and
glutathione pathways suggest that the honeybee is actively attempting to eliminate the
imidacloprid by converting it to a less toxic form. This is crucial, given the reduced
amounts of detoxification genes in honeybees and their unique sensitivity to insecticides
(Appendix Q).
The Longevity Regulating Pathway compares multiple species. For the purpose of
this discussion, the pathway in flies will be considered since they are the most
comparable to the honeybee. In this pathway, Control-4h-Treatment-0h (Appendix P)
comparison shows that various forms of SOD and adenylate cyclase are up-regulated and
heat shock protein cognate (the constitutive form) and insulin-like receptor-like are
down-regulated in imidacloprid-treated bees. The change in regulation of SOD again
supports that a switch occurred between the Treatment-4h and Control-4h; however, there
is no real change in longevity when comparing the treatment group to the control group
Appendix Q). It is interesting to note that HSP70 remained down-regulated (not
statistically significant) while the results of the HSP 70 assay showed a significant
change between concentrations. This could be explained by post-translational regulation,
which is faster than that of post-transcriptional regulation.
The Oxidative Phosphorylation Pathway, located on the cristae of the
mitochondria, consists of a collection of complexes and enzymes responsible for the
aerobic production of ATP. In the imidacloprid treatment group, Complexes I through IV
are up-regulated resulting in an increase in electron transport down the chain and
subsequent increase in the flux of protons across the cristae to the intermembrane space.
79

The last step in the pathway is the ATPase that uses the energy produced by the proton
gradient to synthesize ATP. ATPases are classified into three groups; F, V, and P types.
F-type ATPases are the more common proton-translocating ATP synthases found in the
cristae of mitochondria and also in the membranes of bacteria and chloroplasts that
catalyze the hydrolysis or synthesis of ATP (Ozawa et al. 2000). V-type ATPase use the
energy of hydrolysis to pump protons across membranes to establish electrochemical
gradients that can be used to drive transport processes (Beyenbach and Wieczorek 2006).
The important genes to note are subunit c in both the F-type and V-type ATPase. In the
control group (Appendix P), both these genes are down-regulated (not statistically
significant) while in the imidacloprid-treated group the F-type ATPase is up-regulated
(not statistically significant) while the V-type ATPase is significantly down-regulated
(Appendix Q). Our findings suggest that aerobic metabolism in the honeybee is
significantly impacted at multiple steps in the oxidative phosphorylation pathway;
however, from our findings, whether or not the changes we report increase or decrease
aerobic metabolism cannot be predicted. Other studies suggest that the effects of
neonicotinoid pesticides decrease aerobic metabolism which decreases energy stores
(triglycerides and glycogen) and metabolic fuel (i.e., glucose) (Cook 2019, Tong et al.
2019). Therefore, we can hypothesize that the acute gene changes in oxidative and
glucose metabolism disrupt normal aerobic metabolism and contribute to energetic stress
in honeybees.
A closer look at the function or functional clusters from DAVID revealed
noteworthy categories of enriched clusters of genes that were up-regulated were those
involved in transcription, translation, and pheromone/general odorant binding (Table 7).
Enriched clusters of genes that were down-regulated include transmembrane proteins and
those engaged in intracellular signaling, ATP and GTP binding, ubiquitin-protein
transferase activity, and zinc fingers (Table 7). To determine if a function or functional
cluster was considered to be significant, the FDR had to have a p-value of <0.05. Further
analysis showed an inverse relationship between the enrichment score and FDR p-value;
the higher the enrichment score, the smaller the p-value making that function of the
cluster more enriched. A complete list of unique genes within each annotation cluster can
be found in Appendix R.
80

Table 7: Enriched genes at Treatment-4h and Control-0h. DAVID was used for enrichment analysis.
Genes that had an FDR<0.05 were considered significantly enriched.
Annotation
Cluster
1

2

3

4

5

6

1

2

3
4

6

Cluster Function

Upregulated

Zinc ion binding
Zinc finger, RING/FYVE/PHD-type
Zinc finger, RING-type
RING
Membrane
Transmembrane helix
Transmembrane
Integral component of membrane
Small GTPase mediated signal transduction
Small GTP-binding protein domain
GTP binding
Ubiquitin-protein transferase activity
Ubl conjugation pathway
Ligase helix
HECT
HECTc
Serine/threonine-protein kinase, active site
Protein kinase, ATP binding site
Kinase
Serine/threonine-protein kinase
Protein serine/threonine kinase activity
Protein kinase, catalytic domain
S TKc
Protein kinase-like domain
Nucleotide-binding
ATP-binding
Intracellular signal transduction
Protein kinase C-like, phorbol
ester/diacylglycerol binding
C1
Ribosome
Structural constituent of ribosome
Translation
Ribosome
Sm
Ribonucleoprotein LSM domain
Like-Sm (LSM) domain
Small nucleolar ribonucleoprotein complex
Large ribosomal subunit
Translation protein SH3-like domain
RNA recognition motif domain
Nucleotide-binding, alpha-beta plait
Nucleotide binding
RRM
IG
Immunoglobulin subtype
Immunoglobulin domain
Immunoglobulin-like domain
IGc2
Immunoglobulin subtype 2

81

DownRegulated
√

Enrichment
Score
5.46

√

4.8

√

2.94

√

2.8

√

2.79

√

2.78

√

22.67

√

2.54

√

1.93

√

1.82

√

1.38

7
11
12

RNA polymerase
DNA-directed RNA polymerase activity
Transcription, DNA-templated
Pheromone/general odorant binding protein
PhBP

√

1.27

√
√

0.96
0.95

A comparison of Treatment-4h versus Control-0h and Control-4h versus
Treatment-0h provided a clearer picture of what is happening to pesticide-treated
honeybees at the gene function level. The first noticeable difference is the total amount of
annotation clusters found when the significant genes were inputted into DAVID for
analysis. Control-4h versus Treatment-0h (Appendix S) found only two annotation
clusters in comparison to the Treatment-4h versus Control-0h found multiple annotation
clusters, six of which are relevant to this study (Appendix R); Transmembrane (up- and
down-regulated) and Immunoglobulin (up-regulated).
Another difference between the two groups was that reactive oxygen species, heat
shock protein, and apoptosis were all found to be down-regulated in Treatment-4h versus
Control-0h but the dual oxidase gene (although not considered significant p>0.05) was
down-regulated in the Control-4h versus Treatment-0h. Part of this study looked at SOD
as an indicator of oxidative stress and HSP as an indicator of cellular stress. The dual
oxidase family has the sole function to produce reactive oxidative species (Ameziane-ElHassani et al. 2016). This function correlates with the reduction of the gene in Control-4h
versus Treatment-0h; if the honeybees are not in contact or producing harmful oxygen
molecules then an increase in this gene is unnecessary. As a whole, ROS in the
Treatment-4h versus Control-0h were significantly reduced. In turn resulting to a buildup of toxic oxygen molecules. Although the KEGG analysis showed an elevation in the
pesticide-treated honeybees (Figure 32 and 34) in comparison to the sucrose-treated
honeybees (Appendix P), it was not considered a significant elevation to compensate for
the amount of toxicity within the cells. This finding further provides evidence that both
oxidative and cellular stress are occurring when honeybees do come in contact with
pesticides.
It is apparent that Annotation Cluster 1 in down-regulated genes (Table 7) had the
highest enrichment score. Zinc fingers are proteins that interact with DNA, RNA and
other proteins, and are essential for maintaining homeostasis. Zinc fingers regulate
82

cellular processes such as transcription, ubiquitin-dependent protein degradation, signal
transduction, DNA repair, cell proliferation, differentiation, and apoptosis (Cassandri et
al. 2017). These specific proteins have been implicated in the development of
neurodegenerative disorders. Zinc fingers provide the “checks and balances” within the
honeybee. If these necessary genes are altered for any reason, the “checks and balances”
no longer exist causing a slow downward spiral of control. This study is an excellent
example of that.

VI.

Limitations
Several limitations have been identified in this research. A relatively small sample

of six honeybees per treatment group was sent for transcriptome analysis due to the high
cost of RNA sequencing. However, this number is typical of transcriptome analysis
(Christen et al. 2018) and is in the acceptable sample size recommended by the Keck
Center at the University of Illinois Urbana-Champaign. Due to a switch of the Control-4h
and Treatment-4h samples when labeling, the effects of sublethal dose of imidacloprid on
the individual genes and pathways has to be confirmed by qRT-PCR before publishing
the results of Experiment 2. Despite corroborating evidence for the post-hoc switch
hypothesis, particularly the elevated SOD activity in imidacloprid-treated bees in
Experiment 1, further testing would be required to confirm the results of our
transcriptome analyses.
Another limitation of this research is the type of the post-RNA sequencing
analysis employed. DAVID’s ability to use gene set enrichment to help understand the
meaning behind large lists of genes is valuable; however, there are limitations to what the
database can and cannot do. Of the 10,597 genes, only those genes considered
significantly (p<0.1) up-regulated and down-regulated were uploaded into DAVID
separately. This is because DAVID cannot distinguish the representation by just the
Entrez ID given; the user must tell DAVID what is being uploaded. The pathway analysis
that was used was KEGG, an over-representation analysis (ORA). This method identifies
the pathway membership from the inputted gene list. Some limitations of ORA include:
1) the analysis only considers the number of genes in a pathway and does not account for
over-expression or under; 2) a gene list must be specified and often uses significant hits
83

from the expression analysis, thereby missing more subtle effects; 3) the analysis
assumes each gene is independent of other genes, and assumes each pathway is
independent; and 4) has no statistical output (Farhad et al. 2020). A more inclusive
analysis to use would be Pathway Topology (PT) which uses a functional class scoring
methodology (FCS). FCS calculates gene-level expression statistics from expression data,
uses these statistics to calculate pathway-level statistics, and identifies pathways that are
significantly affected (Bayerlová et al. 2015). Pathway topology incorporates information
about the pathway members such as cellular location and similarity of protein structure
similarity. PT also calculates Perturbation Factor (PF) which is the total effect on a
specific gene and Impact Factor (IF) which is the sum of all PF’s for each gene in a
pathway (Garcίa-Campos et al. 2015). A Pathway Topology program such as
iPathwayGuide could not be utilized due to lack of reference for Apis mellifera.
A final limitation relates to the inherent challenges of conducting toxicology
research on honeybees. By design, this laboratory study evaluated the effects of a single
pesticide on adult bees under controlled conditions. It could be argued that this may not
reflect the exposure to pesticides under field conditions, where bees are subjected to
multiple chemicals, pathogens, varying weather conditions and foraging environments,
and genetic predisposition. There is debate in the literature regarding the imidacloprid
doses used in laboratory studies (Carreck and Rarnieks 2014). Entine (2018) contends
that honeybees in the field have greater opportunities to forage on food sources that are
not contaminated with neonicotinoids, making a case that bee’s dilute pesticides with
“clean forage.” This is in contrast to laboratory bees which are solely fed neonicotinoidspiked “food.”
However, there are several counterarguments to be made. First, clean forage is
difficult to find in agricultural landscapes. For example, neonicotinoid pesticide residues
have been found to contaminate soils and persist in crops and wildflowers in agricultural
landscapes in Europe, where a ban has been in place (Wintermantel et al. 2020). Second,
due to memory and learning, foraging bees display flower constancy and show limited
recruitment to alternate floral resources (Wagner et al. 2013, Van Nest et al. 2018).
Authors have argued that the excellent olfactory sense of bees allows them to detect and
avoid pesticides, but preliminary evidence indicates that this is not true (Plascencia et al.
84

2015). Third, agricultural monocultures provide the largest volume of forage to bees,
providing noticeable hive growth. Such monocultures provide a much higher
contaminated nectar flow than alternate sources can provide of uncontaminated nectar.
Field borders and hedgerows play a role in maintaining persistence of hives when crops
are not in bloom (Sardinas and Kremen 2015, Dolezal et al. 2019). The sublethal doses of
imidacloprid used in this study do represent environmentally relevant exposure rates for
bees foraging in both rural and urban settings. In the wild, even a small loss of motor
coordination, activity level, etc. results in less foraging, increased mortality, less effective
foraging and sublethal stress that can be harmful to the hive (Bryden et al. 2013).
Responses to pesticides may also differ following the oral exposures used in our
study compared to field exposures that may include direct physical contact of bees with
contaminated dust or pollen. Bees in our laboratory study received acute 4h exposures to
imidacloprid. The response of imidacloprid on the transcriptome may differ from chronic
exposures.

VII.

Future Directions

Our research identified 4,205 genes that had experienced significant up- and
down-regulation after acute sublethal exposure to the neonicotinoid pesticide
imidacloprid. To help confirm the imidacloprid-induced changes in the gene expression
that was detected from the RNA-Seq study, future hypothesis-driven research will be
conducted to validate the RNA-seq results. The future study will:
i.

Identify and confirm genes that do not show changes in gene
expression, then utilize the genes as reference genes for gene
expression studies in the honeybee brain tissue.

ii.

Directly test the hypothesis that Treatment-4h and Control-4h
samples were switched and mislabeled prior to shipping to
University of Illinois Urbana-Champaign by confirming gene
expression changes that were detected in our RNA-seq study that
parallel the results from the SOD and HSP70 Assay. The specific
genes that will be tested are the Heat Shock Cognate (HSC70) in

85

the stress apoptotic pathway and Superoxide Dismutase (SOD) in
the oxidative stress pathway.
iii.

Confirm the gene expression changes for other high priority genes
that connect organism-level responses with underlying molecular
mechanisms and biological functions. The key focus will be on:
1. Circadian rhythm gene pathways that correspond with
altered foraging behavior (Karahan et al. 2015, Cakmak et
al. 2018)
2. Cell communication, signal transduction, and drug binding
gene sets that may correspond with motor coordination
(Llewellyn 2020)
3. Response to stimulus and cellular response to stimulus gene
sets that may correspond with performance on sucrose
sensitivity tests (Salazar, In prep.)
4. Macromolecule modification and protein modification gene
sets that may correspond with recovery from or survival of
acute intoxication stress.

VIII. Conclusion
The overall goal of our research is to describe the integrated responses by the
honeybee to sublethal doses of the common neonicotinoid, imidacloprid. The two
objectives of Experiment 1 were: 1) to determine the effects of sublethal doses of
imidacloprid on motor function and overall cellular stress as determine by the levels of
HSP70 and the oxidative stress enzyme SOD, and 2) to identify a sublethal dosage of
imidacloprid that yields a significant cellular stress response to use in transcriptome
studies in Experiment 2. We hypothesized that bees exposed to sublethal doses of the
neonicotinoid imidacloprid will display impaired motor functions and a significant
cellular stress response reflected by increased levels of HSP70 and the oxidative enzyme
SOD. Furthermore, we hypothesized that a dose of imidacloprid could be determined that
yields a significant cellular stress response in honeybee brains. The results of these
experiments confirm numerous literature reports that field-relevant, sublethal doses of
86

imidacloprid impact honeybee physiology and behavior (Pereira et al. 2020, EbanRothschild and Bloch 2011, Honeybee Genome Sequencing Consortium 2006, Wu-Smart
and Spivak 2016, Christen et al. 2018, Dussaubat et al. 2006, Chakrabarti et al. 2014,
Simone-Finstrom et al. 2016). We observed that motor responses of honeybees were
impaired and levels of HSP70 and SOD were elevated when honeybees were acutely
exposed to sublethal doses of imidacloprid. Experiment 1 achieved these goals by
identifying a sublethal dose (0.9 ng/bee, < TD50 for motor impairment, 1.78-2.25 ng/bee)
that coincides with pesticide-induced cellular and oxidative stress, but does not cause
motor impairment. This was important because it is a sublethal dose likely to be
encountered in nature, a dose which bees would appear to be unaffected in an apiary or in
direct-mortality toxicological assays. At the same dose, bees display elevated pesticideinduced, sublethal stress that may contribute to CCD.
Experiment 2 examined the effects of this sublethal dose (0.9 ng/bee) on gene
expression in the honeybee brain. The specific objectives were: 1) to compare overall
gene expression in brain tissue of the control and imidacloprid treatment groups; 2) to
determine which functional classes of genes are up-regulated or down-regulated
following imidacloprid treatment; and 3) to examine difference in gene expression that
occur following imidacloprid treatment in pathways regulating detoxification, heat shock
proteins, oxidative enzymes, energy metabolism, circadian rhythms, cell signaling,
autophagy, apoptosis, and longevity. We hypothesized that gene expression patterns will
be altered in the brain tissue in imidacloprid treated bees compared to controls and will
express different functional classes of genes compared to control bees. We proposed that
pathways regulating detoxification, heat shock proteins, oxidative enzymes, energy
metabolism, and apoptosis would be up-regulated, while pathways regulating circadian
rhythms, cell signaling, and longevity would be down-regulated in imidacloprid-treated
bees.
A total of 4,205 genes were identified as differentially expressed. However,
differential expression analysis, a multi-dimensional scaling plot, and the one-way
ANOVA heat map revealed that the greatest changes in gene expression occurred in the
Control-4h group, contrary to our hypothesis and previous transcriptome studies in the
literature. Our independent results in Experiment 1 and 2 provide corroborating evidence
87

supporting the post-hoc hypothesis that the Control-4h and the Treatment-4h samples
were switched when the samples were mislabeled prior to shipping to the University of
Illinois Urbana-Champaign. Future studies are planned to directly test the post-hoc
hypothesis by confirming gene expression changes by RT-qPCR that were detected in our
RNA-sequencing study that parallel the results from the SOD and HSP70 Assay. The
specific genes that will be tested are the Heat Shock Cognate 70 (HSC70) in the stress
apoptotic pathway and Superoxide Dismutase (SOD) in the oxidative stress pathway.
Although a switch of the samples likely occurred, the comparison between the
Control-4h and the Treatment-4h groups remain valid, since only the direction of
differential gene expression (up-regulated versus down-regulated) would be affected.
Pathway analysis indicated that multiple pathways were affected including oxidative
phosphorylation, longevity regulating pathway, apoptosis, peroxisome, FOXO signaling,
drug metabolism-cytochrome P450, metabolism of xenobiotics by cytochrome P450,
circadian rhythm, and glutathione metabolism. Similar pathways have been altered in
bees exposed to nicotine (du Rand et al. 2015), ethanol (Hranitz et al. unpublished), and
thiamethoxam (Shi et al. 2017). The specific alterations we observed support our
prediction that pathways that regulate detoxification, heat shock proteins, oxidative
enzymes, energy metabolism, and apoptosis will be up-regulated while pathways that
regulate circadian rhythms, cell signaling, and longevity will be down-regulated in
imidacloprid-treated bees. These alterations in gene networks relate to key biological
functions of honeybees that have the potential to affect the viability of the colony.
Long term, our aim is to understand the specific gene networks and cellular
pathways affected by sublethal imidacloprid intoxication in order to understand the
molecular mechanisms of neonicotinoid toxicity and its contribution to pollinator
declines. This research serves as a foundation for future hypothesis-driven gene
expression studies that relate specific molecular changes that occur with neonicotinoid
toxicity to biological functions. Future studies can target gene expression changes for
other high priority genes that connect organism-level responses in the field with
underlying molecular mechanisms. This will not only allow us to understand
neonicotinoid toxicity but may also prevent pollinator declines and CCD.

88

IX.

Literature Cited

Abramson CI, Kandolf A, Sheridan A, Donohue D. Bozic J, Meyers JE, Benbassat D.
(2004a). Development of ethanol model using social insects: III. Preferences for
ethanol solutions. Physiological reports 94(1):227.
Abramson CI, Place AJ, Aquino IS, Fernandez A. (2004b). Development of ethanol
model using social insects: IV. Influence of ethanol on the aggression of Africanized
Honeybees (Apis mellifera L.) Physiological Reports 94(3): 1107.
Abramson CI, Stone SM, Ortez RA, Luccardi A, Vann KL, Hanig KD, Rice J. (2000).
Development of an ethanol model using social insects 1: Behavior studies of the
Honeybee (Apis mellifera L.). Alcoholism-Clinical and Experimental Research
24(8):1153.
Abramson CI, Fellows FW, Browne BL, Lawson A, Ortiz RA. (2003). Development of
an ethanol model using social insects: II. Effects of Antabuse ® on consumatory
responses and learned behavior of the Honeybee (Apis mellifera L.). Physiological
Reports 92(2): 365.
Abramson CI, Sanderson C, Painter J, Barnett S, Wells H. (2005). Development of an
ethanol model using social insects: V. Honeybee foraging decisions under the influence
of alcohol. Alcohol 36(3):187.
Agriculture and Consumer Protection. (1990). Chapter 2 Colony life and social
organization. www.fao.org.
Albaugh LLC. (2017). Company. Retrieved from albaughllc.com/company/.
Ameziane-El-Hassani R, Schlumberger M, Dupuy C. (2016). NADPH oxidases: new
actors in thyroid cancer? Nat Rev Endocrinol. 12(8):485-94. doi:
10.1038/nrendo.2016.64.
Atmowidjojo, A. H., Wheeler, D. E., Erickson, E. H., & Cohen, A. C. (1997).
Temperature tolerance and water balance in feral and domestic Honeybees, Apis
mellifera L. Comparative Biochemistry and Physiology Part A: Physiology, 118(4),
1399-1403.
Atwood, D and Paisley-Jones, C. (2017). Pesticide Industry Sales and Usage 2008-2012,
Market Estimate. epa.gov
Alavanja M. C. (2009). Introduction: pesticides use and exposure extensive
worldwide. Reviews on environmental health, 24(4), 303–309.
https://doi.org/10.1515/reveh.2009.24.4.303

89

Barthell JF, Hranitz JM, Thorp RW, and Shue MK. (2002). High temperature
responses in two exotic leafcutting bee species: Megachile apicalis and M. rotundata
(Hymenoptera: Megachilidae). The Pan-Pacific Entomologist. 78(4): 235-246.
Bayer Crop Science. (2017). Insecticides. Retrieved from
https://www.cropscience.bayer.us/products/insecticides/.
Bayerlová, M., Jung, K., Kramer, F. et al. (2015). Comparative study on gene set and
pathway topology-based enrichment methods. BMC Bioinformatics 16: 334.
https://doi.org/10.1186/s12859-015-0751-5
Bee Culture. (2019). Catch the Buzz- during the 2018-2019 winter an estimated
37.7% of managed Honeybee colonies in the United States were lost.

https://www.beeculture.com/catch-the-buzz-during-the-2018-2019-winter-an-estimated-37-7of-managed-honey-bee-colonies-in-the-united-states-were-lost/

Benbrook, C. (2009). Impacts of genetically engineered crops on pesticide use: the
first thirteen years. Boulder, CO. The Organic Center.
Berlett BS and Stadtman ER. (1997). Protein oxidation in aging, disease, and oxidative
stress. J Biol Chem 272(33):20313-6.
Beyenbach, Klaus and Wieczorek, Helmut. (2006). The V-type H+ ATPase: molecular
structure and function, physiological roles and regulation. Journal of Experimental
Biology 209:577-589.
Blacquière T, Smagghe G, van Gestel C, Mommaerts V. (2012). Neonicotinoids in bees:
a review on concentrations, side-effects and risk assessment. Ectotoxicology 21:973992. DOI 10.1007
Bozic J, Abramson CI, Bedencic M. (2006). Reduced ability of ethanol drinkers for
social communication in Honeybees. (Apis mellifera carnica Poll.). Alcohol 38(3):179.
Bozic J, DiCesare J, Wells H, Abramson CI. (2007). Ethanol levels in Honeybee
hemolyph resulting from alcohol ingestion. Alcohol 41(4):281.
Bryden J., Gill RJ, Mitton RAA, Raine NE, Jansen VAA. (2013). Chronic sublethal
stress causes bee colony failure. Ecology Letters 16(12):1463-1469.
Cabirol, A., & Haase, A. (2019). The Neurophysiological Bases of the Impact of
Neonicotinoid Pesticides on the Behaviour of Honeybees. Insects, 10(10), 344.
https://doi.org/10.3390/insects10100344
Cakmak I, Hranitz J, Blatzheim L, Bower C, Polk T, Levinson B, Wells H. (2018).
Effects of thiamethoxam on the behavior of foraging Honeybees with artificial
flower choices. Uludag Bee Journal 18(1):2-13.
90

Calabrese EJ, Stanek EJ, III, Nascarella MA, Hoffman GR. (2008). Hormesis Predicts
Low-Dose Responses Better Than Threshold Models. International Journal of
Toxicology 27(5):369.
Carreck, Norman and Ratnieks, Francis. (2014). The dose makes the poison: have
“field realistic” rates of exposure of bees to neonicotinoid insecticides been
overestimated in laboratory studies?. Journal of Apicultural Research. 53(5):607-614.
DOI: 10.3896/IBRA1.53.5.08
Cassandri M, Smirnov A, Novelli F, Pitolli C, Agostini M, Malewicz M, Melino G, and
Rashellà G. (2017). Zinc-finger proteins in health and disease. Cell Death Discovery.
3.
Chacon-Almeida L, Simões Z, and Bitondi M. (2000). Induction of heat shock proteins
in the larval fat body of Apis mellifera L. bees. Apidologie, Springer Verlag. 31(4):
487-501
Chakrabarti P, Rana S, Sarkar S, Smith B, Basu P. (2015). Pesticide-induced oxidative
stress in laboratory and field populations of native Honeybees along intensive
agricultural landscapes in two eastern indian states. Apidologie 46(1):107-29.
Christen V, Schirrmann M, Frey J, and Fent K. (2018). Global Transcriptomic Effects
of Environmentally Relevant Concentrations of the Neonicotinoids Clothianidin,
Imidacloprid, and Thiamethoxam in the Brain of Honeybees (Apis mellifera).
Environmental Science and Technology. 2(13): 7534-7544.
Claudianos, C., Ranson, H., Johnson, R. M., Biswas, S., Schuler, M. A., Berenbaum, M.
R., Feyereisen, R., & Oakeshott, J. G. (2006). A deficit of detoxification enzymes:
pesticide sensitivity and environmental response in the honeybee. Insect molecular
biology, 15(5), 615-636. https://doi.org/10.1111/j.1365-2583.2006.00672.x
Cook SC. (2019). Compound and Dose-Dependent Effects of Two Neonicotinoid
Pesticides on Honey Bee (Apis mellifera) Metabolic Physiology. Insects, 10(1), 18.
https://doi.org/10.3390/insects10010018
Coulon, M., Schurr, F, Martel, AC, Cougoule, N, Begaud, A., Mangoni, P, Dalmon, A,
Alaux, C. LeConte, Y, Thiery, R, Ribiere-chabert, M, Dubois, E. Metabolization of
thiamethoxam (a neonicotinoid pesticide and interaction with the chronic bee
paralysis virus in honeybees. Pesticide Biochemistry and Physiology, Vol. 144, p10-18,
2018.

91

Decourtye A, Lacassie E, Pham-Delegue MH. (2003). Learning performances of
honeybees (Apis mellifera L.) are differentially affected by imidacloprid according to
the season. Pest Manag Sci 59:269-278.
Decourtye A, Devillers J, Cluzeau S, Charreton M, Pham-Delegue MH. (2004a). Effects
of imidacloprid and deltamethrin on associative learning in honeybees under semifield and laboratory conditions. Ecotoxicol Environ Saf 57:410-419.
Decourtye A, Armengaud C, Renou M, Devillers J, Cluzeau S, Gautheir M, PhamDelegue MH. (2004b). Imidacloprid impairs memory and brain metabolism in the
honeybee (Apis mellifera L.). Pest Biochem Physiol 78:83-92.
Deng C, Graham R, and Shukla R. (2001). Detecting and Estimating Hormesis Using a
Model-Based Approach. Human and Ecological Risk Assessment: An International
Journal. 7: 849-866.
DiPrisco G, Cavaliere V, Annoscia D., Varricchio P, Caprio E, Nazzi,F, Gargiuli
G,Pennacchio F. (2013). Neonicotinoid clothianidin adversely affects insect immunity
and promotes replication of a viral pathogen in Honeybees. Proc of the Natl Acad Sci.
110(46):18466-18471.
Dolezal AG, St Clair AL, Zhang G, Toth AL, and O’Neal ME. (2019). Native habitat
mitigates feast-famine conditions faced by Honeybees in an agricultural landscape.
Proceedings of the National Academy of Sciences of the United States of America.
116(50): 25147-25155.
Drake, J. M. & Kramer, A. M. (2011). Allee Effects. Nature Education
Knowledge 3(10):2.
du Rand EE, Smit S, Beukes M, Apostolides Z, Pirk C, Nicolson S. (2015).
Detoxification mechanisms of Honeybees (Apis mellifera) resulting in tolerance of
dietary nicotine. Scientific Reports 5:11779 DOI: 10.1038.
Dussaubat C, Maisonnasse A, Crauser D, Tchamitchian S, Bonnet M, Cousin M,
Kretzschmar A, Brunet J, Le Conte Y. (2016). Combined neonicotinoid pesticide and
parasite stress alter honeybee queens’ physiology and survival. Scientific Reports
6:31430.
Eban-Rothschild A., Bloch G. (2012). Circadian Rhythms and Sleep in Honeybees.
In: Galizia C., Eisenhardt D., Giurfa M. (eds) Honeybee Neurobiology and Behavior.
Springer, Dordrecht
Ellis JD and Mortenson AN. (2011). Featured Creatures. University of Florida:
Entomology and Nematology.
http://entnemdept.ufl.edu/creatures/misc/bees/cape_honey_bee.htm
92

Engel, MS. (1999). The taxonomy of recent and fossil Honeybees (Hymenoptera:
Apidae: Apis). Journal of Hymenoptera Research. 8:165-195.
Entine, Jon. (2018). “Gold standard” assessing neonicotinoids: Field bee hive studies
find pesticides not major source of health issues. Genetic Literacy Project.
Farhad M, Ovens K, Hogan DJ, and Kusalik AJ. (2020). Gene Set Analysis: Challenges,
Opportunities, and Future Research. Frontiers in Genetics. 11
Farhan, M., Wang, H., Gaur, U., Little, P. J., Xu, J., & Zheng, W. (2017). FOXO
Signaling Pathways as Therapeutic Targets in Cancer. International journal of
biological sciences, 13(7), 815–827.
Feder ME, Blair N, Figeuras H. (1997). Natural thermal stress and heat-shock protein
expression in Drosophila larvae and pupae. Functional Ecology 11:90-100.
Feder ME and Hoffman G. (1999). Heat-shock proteins, molecular chaperones, and
the stress response: evolutionary and ecological physiology. Annual Review of
Physiology. 61(1): 24-282.
Finkel T and Holbrook NJ. (2000). Oxidants, oxidative stress and the biology of
ageing. Nature 408(6809):239.
Frazier M, Mullin, C., Frazier, J. Ashcraft, S. (2008) What have pesticides got to do
with it? American Bee Journal. Vol 148, p521-523.
Garcίa-Campos M, Espinal-Enrίquez J, and Hernández-Lemus E. (2015). Pathway
Analysis: State of the Art. Front. Physiol. DOI:10.3389.
Gervais, J. A., Luukinen, B., Buhl, K., & Stone, D. (2010). Imidacloprid general fact
sheet. Oregon State University Extension Services: National Pesticide Information
Center.
Girolami, V., Mazzon, L., Squartini, A., Mori, N., Marzaro, M., Di Bernardo, A., et al.
(2009). Translocation of neonicotinoid insecticides from coated seeds to seedling
guttation drops: A novel way of intoxication for bees. Journal of Economic
Entomology, 102(5), 1808-1815.
Goulson, David. (2013). An overview of the environmental risks posed by
neonicotinoid insecticides. Journal of Applied Ecology 50, 977-987.
Henry M, Reguin M, Requier F, Rollin O,Odoux J-F, Aupinel P, Qiao D, Seidler F,
Slotkin T. (2005). Oxidative mechanisms contributing to the developmental
neurotoxicity of nicotine and chlorpyrifos. Toxicol and Applied Pharmacol. 206(1):1726.
93

Honeybee Genome Sequencing Consortium. (2006). Insights into social insects from
the genome of the honeybee Apis mellifera. Nature 443(26): 931-949.
Honeymell. (2015). Reproduction of Honeybees.
https://www.honeymell.com/en/bees/25-reproduction-of-honey-bees

Hrantitz et al. (2010). An ethanol-induced hermetic stress response in Honeybees
(Apis mellifera) brain tissue. Integrative and Comparative Biology 50:E245-E245.
Hranitz JM and Barthell JF. (2003). Heat shock protein 70 during development of the
leafcutting bee, Megachile rotundata (Hymenoptera: Megachilidae). Southwestern
Entomologist. 28(2): 97-104.
Hranitz JM, Barthell JF, Abramson CI, Brubaker KD, Wells H. (2009a). Stress protein
responses in Honeybees: Is it useful to measure stress responses of individual bees in
the hive? Uladag Bee Journal 9(2): 60-71.
Hranitz JM, Abramson CI, Carter RP. (2009). Ethanol increases HSP70 concentrations
in honeybee (Apis mellifera L.) brain tissue. Alcohol 44:275-282.
Hunt, GJ. (2007). Flight and fight: A comparative view of the neurophysiology and
genetics of Honeybee defensive behavior. Journal of Insect Physiology 53(5):399-410.
Ighodaro OM and Akinloye OA. (2018). First line defense antioxidants-superoxide
dismutase (SOD), catalase (CAT) and glutathione peroxidase (GPX): Their
fundamental role in the entire antioxidant defense grid. Alexandria Journal of
Medicine 54(4):287-293.
Illumina. (2010). Illumina Sequencing Technology.
https://www.illumina.com/documents/products/techspotlights/techspotlight_sequencing.p
df
Johnson, B.R. (2008). Within-nest temporal polyethism in the Honeybee. Behav Ecol
Sociobiol 62:777.
Johnson RM, Wen Z, Schuyker MA, Berenbaum MR. (2006). Mediation of pyrethroid
insecticide toxicity to hyoney bees (Hymenoperia:Apidae) by Cytochrome P450
monooxygenases. Econ Entomol, Aug 99 (4) 1046-50
Jones JC, Helliwell P, Beekman M, Maleszka R, and Oldroyd BP. (2005). The effects of
rearing temperature on developmental stability and learning and memory in the
Honeybee, Apis mellifera. Journal of Comparative Physiology a-Neuroethology Sensory
Neural and Behavioral Physiology 191(12):1121-1129.
Karahan A, Cakmak I, Hranitz JM, Karaca I, Wells H. (2015). Sublethal imidacloprid
effects on Honeybee flower choices when foraging. Ecotoxicology 24:2017-2025.
94

KEGG Mapper. (2020). https://www.genome.jp/kegg/mapper.html
Kessler, S., Tiedeken, E., Simcock, K. et al. (2015). Bees prefer foods containing
neonicotinoid pesticides. Nature 521, 74–76.
Kuropatnicki Andrzej, Szliszka Ewelina, and Krol Wojciech. (2013). Historical Aspects
of Propolis Research in Modern Times. Evid Based Complement Alternat Med. DOI:
10.1155/2013/964149.
Le Conte Y, Ellis M, and Ritter W. (2010). Varroa mites and Honeybee health: can
Varroa explain part of the colony losses? Apidologie 41(3):353-363. DOI:
https://doi.org/10.1051/apido/2010017.
Maiese, Kenneth. (2015). FOXO Proteins in the Nervous System Analytical Cellular
Pathology. 2015. DOI:10.1155/2015/569392
Mao W, Schuler M, and Berendaum M. (2013). Honey constituents up-regulate
detoxification and immunity genes in the western Honeybee Apis mellifera. PNAS
110(22): 8842-8846.
Mater Methods. (2012). 2:201.
Merritt, Chris. (2011-2020). Hormesis: Why you may not be reaching your goals?
https://www.bspnova.com/hormesis-why-you-might-not-be-reaching-your-goals/
Möller, Gregory. (NA). Dose-Response Relationships. University of Idaho.
https://www.webpages.uidaho.edu/foodtox/lectures/lecture05/L5-Dose%20%20Response%20Relationships.pdf
National Center for Biotechnology Information. (accessed on July 22,2020). PubChem
Database. CID=86287518,
https://pubchem.ncbi.nlm.nih.gov/compound/Imidacloprid
Nature Protocols. (2009). 4(1):44. https://david.ncifcrf.gov/
Ozawa, K., Meikari, T., Motohashi, K., Yoshida, M., & Akutsu, H. (2000). Evidence for
the presence of an F-type ATP synthase involved in sulfate respiration in
Desulfovibrio vulgaris. Journal of bacteriology, 182(8), 2200–2206.
https://doi.org/10.1128/jb.182.8.2200-2206.2000
Pereira NC, Diniz TO, and Ruvolo-Takasusuki MCC. (2020). Sublethal effects of
neonicotinoids in bees: a review. Sci. Elec. Arch 13(7). DOI:
http://dx.doi.org/10.36560/13720201120.
Pilatic, Heather. (2012). Pesticides and Honeybees: State of Science. www.panna.org

95

Plascencia M, Carson R, Hranitz JM, Barthell JF, Cakmak I, Gonzalez VH. (2015).
Testing the pesticide avoidance hypothesis by bees and flies in a Mediterranean
arthropod community. Integrative and Comparative Biology. 55:E315-E315.
Poquet Y, Vidau C, Alaux C. (2015). Modulation of pesticide response in Honeybees.
Apidologie 47:412-416.
Qiagen. (2020). RNeasy Micro Handbook. HB 1920-002.
Sahebzadeh, Najmeh & Lau, Wei. (2017). Expression of heat-shock protein genes in
Apis mellifera meda (Hymenoptera: Apidae) after exposure to monoterpenoids and
infestation by Varroa destructor mites (Acari: Varroidae). European Journal of
Entomology. 114. 195-202. 10.14411/eje.2017.024.
Salazar T, Young C, Naranjo S, Pastor M, Gunes N, Cakmak I, Barthell J, Hranitz J. (In
Prep.) Effects of sublethal doses of insecticides on Honeybees (Hymenoptera:
Apidae). Ecotoxicology.
Sardinas HS and Kremen C. (2015). Pollination services from field-scale agricultural
diversification may be context-dependent. Agriculture Ecosystems and Environment.
207:17-25.
Schroeder ME and Flattum, RF. (1984). The mode of action and neurotoxic properties
of the nitromethylene heterocycle insecticides. Pesticide Biochemistry and Physiology
(22)2:148-160
Seeley, TD. (1985). Honeybee Ecology. Princeton University Press, Princeton, NJ.
Sen Sarma M, Rodriguez-Zas SL, Hong F, Zhong S, Robinson GE. (2009).
Transcriptomic profiling of central nervous system regions in three species of
Honeybee during dance communication behavior. PLoS ONE 4(7):e6408.
Simon-Delso, N., Amaral-Rogers, V., Belzunces, L. P., Bonmatin, J., Chagnon, M.,
Downs, C., et al. (2015). Systemic insecticides (neonicotinoids and fipronil): Trends,
uses, mode of action and metabolites. Environmental Science and Pollution Research,
22(1), 5-34.
Simone-Finstrom M, Li-Byarlay H, Huang MH, Strand MK, Rueppell O, Tarpy DR.
(2016). Migratory management and environmental conditions affect lifespan and
oxidative stress in Honeybees. Scientific Reports 6:32023.
Shi, TF, Wang YF, Liu F, Qi, L, Yu, LS. (2017). Influence of the neonicotinoid
insecticide thiamethoxam on miRNA Expression in the Honeybee (Hymenoptera:
Apidae). Journal of Economic Entomology. Col xx, p 1-7

96

Stoner KA and Eitzer EB. (2012). Movement of soil-applied imidacloprid and
thiamethoxam into nectar and pollen of squash (Cucurbita pepo). PLoS One 7(6)
Sylvester, H., Rinderer, T., & Bolten, A. (1983). Honey sac contents: A technique for
collection and measurement in foraging Honeybees (hymenoptera: Apidae). Journal
of Economic Entomology, 76(1), 204-206.
Tan, K., Chen, W., Dong, S., Liu, X., Wang, Y., & Nieh, J. C. (2014). Imidacloprid
alters foraging and decreases bee avoidance of predators. PLoS One, 9(7), e102725.
Tasman K, Rands S, and Hodge J. (2020). The neonicotinoid insecticide imidacloprid
disrupts bumblebee foraging rhythms and sleep. bioRxiv.
https://www.biorxiv.org/content/10.1101/2020.04.07.030023v1.full
Technical Note: Sequencing: Quality Scores for Next-Generation Sequencing. Assessing
sequencing using Phred quality scoring. (2011). Illumina Inc.
https://www.illumina.com/documents/products/technotes/technote_Q-Scores.pdf
Toma, D. P., Bloch, G., Moore, D., & Robinson, G. E. (2000). Changes in period
mRNA levels in the brain and division of labor in Honeybee colonies. Proceedings of
the National Academy of Sciences of the United States of America, 97(12), 6914–6919.
https://doi.org/10.1073/pnas.97.12.6914
Tomizawa M, Lee DL, and Caasida J. (2000). Neonicotinoid Insecticides: Molecular
Features Conferring Selectivity for Insect versus Mammalian Nicotinic Receptors.
Journal of Agricultural and Food Chemistry 48(12):601624. DOI: 10.1021/jf000873c.51.48
Tong L, Nieh JC, and Tosi S. (2019). Combined nutritional stress and a new systemic
pesticide (flupyradifurone, Sivanto®) reduce bee survival, food consumption, flight
success, and thermoregulation. Chemosphere, 237, 124408.
https://doi.org/10.1016/j.chemosphere.2019.124408
United States Environmental Protection Agency. (2016). Colony Collapse Disorder.
www.epa.gov.
U.S. National Library of Medicine. (2015). PTEN Gene.
https://ghr.nlm.nih.gov/gene/PTEN#:~:text=Normal%20Function,or%20in%20an%20unc
ontrolled%20way.
Valizadegan N and Drnevich J. (2019). John M Hranitz RNA-Seq Report. University
of Illinois.
van Dooremalen, C, Cornelissen, B, Poleij-hok-ahin, C, Blacquiere, T. (2018). Single
and interactive effects of Varroa destructor, Nosema spp., and imidacloprid on
Honeybee colonies (Apis mellifera). Exosphere. Vol 9. pe02378.
97

vanEngelsdorp D, Evans JD, Saegerman C, Mullin C, Haubruge E, Nguyen BK et al.
(2009). Colony collapse disorder: A descriptive study. PLoS ONE 4(8): e6481.
Van Nest BN, Otto MW, and Moore D. (2018). High experience levels delay
recruitment but promote simultaneous time-memories in Honeybee foragers.
Journal of Experimental Biology. 221(23).
Van Rossum, G and Drake Jr, FL. (1995). Python reference manual. Centrum voor
Wiskunde en Informatica Amsterdam.
Vannette R, Mohamed A, Johnson B. (2015). Forager bees (Apis mellifera) highly
express immune and detoxification genes in tissues associated with nectar
processing. Scientific Reports 5:16224 DOI 10.1038
Wagner AE, Van Nest BN, Hobbs CN, Moore D. (2013). Persistence, reticence and the
management of multiple time memories by forager Honeybees. Journal of
Experimental Biology. 216(7): 1131-1141
Wainselboim AJ et al. (2000). Trophallaxis in the honeybee Apis mellifera (L.): the
interaction between flow of solution and sucrose concentration of the exploited food
sources. Anim Behav 59(6):1177-1185.
Wilson, EO. (1971). The Insect Societies. Cambridge, MA: Belknap Press of Harvard
University Press.
Wintermantel D, Odoux JF, Decourtye A, Henry M, Allier F, and Bretagnolle V. (2020).
Neonicotinoid-induced mortality risk for bees foraging on oilseed rape nectar
persists despite EU moratorium. Science of the Total Environment. 704.
Wu MC, Chang YW, Lu KH, and Yang EC. (2012). Gene expression changes in
Honeybees induced by sublethal imidacloprid exposure during the larval stage.
Insect biochemistry and molecular biology 88.
Wu-Smart, J and Spivak, M. (2016). Sub-lethal effects of dietary neonicotinoid
insecticide exposure on Honeybee queen fecundity and colony development Scientific
Reports 6:32108 DOI:10.1038.
Yan H, Jia H, Gao H, Guo X, and Xu B. (2013). Identification, genomic organization,
and oxidative stress response of a sigma class glutathione S-transferase gene
(AccGSTS1) in the Honeybee, Apis cerana cerana. Cell Stress Chaperones. 2013 Jul;
18(4): 415–426. DOI: 10.1007/s12192-012-0394-7
Zemene M, Bogale B, and Derso S. (2015). A Review on Varroa Mites of Honeybees.

98

Zhiguo L, Tiantian Y, Chen Y, Matthew H, Jingfang H, Jingnan HU, Hongyi N, Songkun
S. (2019). Brain transcriptome of Honeybees (Apis mellifera) exhibiting impaired
olfactory learning induced by a sublethal dose of imidacloprid. Pesticide
Biochemistry and Physiology. 156: 36-43.

99

Appendix A: Animal Research (IACUC) Approval.

384784

100

Appendix B: Literature data on neonicotinoid residues in bee-collected pollen, honey, and bees (Blacquière et al. 2012).
Substrate

Country

Dominant
Crop

Design and
Analysis

Chemical

Sample
Size

%
Positive

Concentration
(ug kg-1 fresh
weight)

Notes

References

Many
pesticides
detected

Mullin et al.
(2010)

<2.0
<1.0
<2.0
<1.0

Many
pesticides
detected

Mullin et al.
(2010)

- |2.4-13.6
<1.0
<5.0
0.1|<1.0-8.0

Many
pesticides
detected

Mullen et al.
(2010)

Mean | +
Range
USA(Florida,
California,Pennsylvania
or 13 states) + samples
from outside USA and
Canada

No data

Honey
Bees

USA(Florida,
California,Pennsylvania
or 13 states) + samples
from outside USA and
Canada

No data

Bee Wax

USA (Florida,
California, Pennsylvania
or 13 states) + samples
from outside USA and
Canada

No data

101

Pollen

2007/2008 survey
included 13 apiaries
in
Florida+California,
47 colonies in
Pennsylvania
orchards and from
“other” samples; no
further data analysis
2007/2008 survey
included 13 apiaries
in
Florida+California,
47 colonies in
Pennsylvania
orchards and from
“other” samples; no
further data analysis
2007/2008 survey
included 13 apiaries
in
Florida+California,
47 colonies in
Pennsylvania
orchards and from
“other” samples; no
further data analysis

Imida.
Thiam.
Aceta.
Thiac.

350

2.9
0.29
3.1
5.4

Imida.
Thiam.
Aceta.
Thiac.

140

0
0
0
0

Imida.
Thiam.
Aceta.
Thiac.

208

0.96
0
0
1.9

101

3.1 | <2.0-912
53.3
1.9 | <5.0-124
5.4 | <1.0-115

Appendix C: Lethal and sub-lethal effects by imidacloprid to individual (organism level) honey bees as determined in different studies by oral
exposure under laboratory conditions (Blacquière et al. 2012).
Species
Apis mellifera
Apis mellifera
Apis mellifera

Exposure
Lab+Oral: acute exposure to 0.12
and 12 ng/ bees
Lab+Oral: acute exposure to 0.23.2 mg/L
Lab+Oral: chronic exposure (no
information on concentration)

102
Apis mellifera

Apis mellifera
Apis mellifera

Lab+Oral: acute exposure to 0.181 ng/bee

Lab+Oral: acute exposure to 0.7
mg/seed
Lab+Oral: chronic exposure (39
days) to sunflower nectar
contaminated with 0.002-0.02
ug/kg-1

Side-effects
Reduction of associative learning
at 12 ng/bee-1
LD50 -48h: 30 ng/bee-=1
Lowest observed effect
concentration on survival of
winter bees: 24 ug/kg-1
Lowest observed effect
concentration on associative
learning via proboscis extension
reflex assay on winter bees (12
ng/bee-1) and summer bees (12
ng/bee-1)
LD50-48h: between 41 and >81
ng/bee-1
No-observed effect dose :</= 1.25
ng/bee; reduced sucrose uptake
by 33% at 81 ng/bee-1
LD50-48h: 4-41 ng/bee-1
No-observed effect concentration
for mortality, feeding activity,
wax comb production, breeding
performance and colony vitality:
0.02 ug/kg-1

102

References
Decourtye et al. (2004a, b)
Decourtye et al. (2003)
Decourtye et al. (2003)

Nauen et al. (2001)

Schmuck et al. (2001)
Schmuck et al. (2001)

Appendix D: Concentrations of imidacloprid causing lethal and sublethal effects on (micro)-colony level in honey bees as determined in different
studies by oral exposure under laboratory conditions (Blacquière et al. 2012).
Species
Apis mellifera

Exposure
Lab+Oral: 100-1,000 ug/L-1

Apis mellifera

Lab+Oral: 0.12-12 ng/bee-1 in
syrup (acute)
Lab+Oral: 24 ug/kg-1 in syrup
(24h)

103
Apis mellifera

Lab+Oral: >100ug/kg-1 (chronic)

Toxicity
500-1,000 ug/L-1: bees
disappeared at the hive/feeding
site up to 24 hours
Increase of the cytochrome
oxidase labeling, negative effect
on the proboscis extension reflex
assay with 12 ng/bee-1 but not
with 0.12 ng/bee-1
Negative effect on the proboscis
extension reflex assay
No-observed effect concentration
survival: 2-20 ug/kg-1 in
sunflower nectar
20 ug/L-1: decrease in foraging
activity; >100 ug/L-1: reduce in
foraging behavior for 30-60
minutes
>50 ug/L-1: increase in interval
between successive visits at a
feeder

103

References
Bortolottie et al. (2003)

Decourtye et al. (2004a, b)

Schmuck (1999); Schmuck et al.
(2001)

Appendix E: Macho® Stock solution preparation.
Stock Solution Preparation
Sucrose Solution:
Molar Mass: 342.3 g/mol
C12H22O11

104

1.5 𝑀𝑀 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 × 342.3

𝑔𝑔
𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 = 513.45 𝑔𝑔 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 × 0.1 = 51.345 𝑔𝑔 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 𝑖𝑖𝑖𝑖 100 𝑚𝑚𝑚𝑚 𝑈𝑈𝑈𝑈𝑈𝑈𝑈𝑈𝑈𝑈 𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 𝑊𝑊𝑊𝑊𝑊𝑊𝑊𝑊𝑊𝑊
𝑚𝑚𝑚𝑚𝑚𝑚

This solution was used as part of the stock solutions for imidacloprid treatment and was refrigerated at 4°C.
Macho Solution:
18,14.37𝑔𝑔
0.4793𝑔𝑔 109 𝑛𝑛𝑛𝑛
𝑚𝑚𝑚𝑚
𝑛𝑛𝑛𝑛
4𝑙𝑙𝑙𝑙𝑙𝑙
=
=
�
� � 3 � = 479,306
𝑗𝑗𝑗𝑗𝑗𝑗
1𝑔𝑔𝑔𝑔𝑔𝑔 3785.41𝑚𝑚𝑚𝑚
𝑚𝑚𝑚𝑚
𝑔𝑔
10 𝜇𝜇𝜇𝜇
𝜇𝜇𝜇𝜇
𝐶𝐶1𝑉𝑉1 = 𝐶𝐶2𝑉𝑉2

𝑛𝑛𝑛𝑛
𝑛𝑛𝑛𝑛
� 𝑉𝑉1 = �3.6
� (10 𝑚𝑚𝑚𝑚)
𝜇𝜇𝜇𝜇
10𝜇𝜇𝜇𝜇
𝑉𝑉1 = 0.0075 𝑚𝑚𝑚𝑚 𝑓𝑓𝑓𝑓𝑓𝑓 25 𝑚𝑚𝑚𝑚 = 0.0187 𝑚𝑚𝑚𝑚
�479,306

Stock 1/5 LD50:
100 μL from Macho jug added to 1,000 mL Ultra Pure Water = 187 μL of stock added to 25 mL stock Sucrose Solution
*Macho is kept at room temperature per pesticide directions.
104

Appendix F: Treatment group preparation.
Treatment
(ng/bee)
1/5 LD50
1/10 LD50
1/20 LD50
1/50 LD50
1/100 LD50
1/500 LD50
Total

[ng/μL] for 10 μL
Feeding
3.6
1.8
0.9
0.36
0.18
0.036

1/5 Stock Solution
(mL)
10
5
2.5
1
0.5
0.005
19.005 mL

All treatment solutions should be stored in the refrigerator at 4℃.

105

1.5 M Sucrose
Solution
0
5
7.5
9
9.5
9.995
40.995 mL

Appendix G: Bovine stock albumin standard preparation.

The BSA standards, determined using C1V1=C2V2 with a stock solution of 7.1 μg/μL of
Bovine Stock Albumin, used are seen in the table below:

Standards

[Protein]

BSA (μL)

Homogenization
Buffer (μL)

1

0.22187

6.25

193.75

2

0.44375

12.5

187.5

3

0.8875

25

175

4

1.775

50

150

5

3.55

100

100

6

5.0

140.8

59.2

7

7.1

200

0

106

Appendix H: Stock bovine HSP70 standard preparation.

Standard

1:10 Standard

[Final] ng

Volume

Volume

HSc70 (μL)

Buffer (μL)

10

4

196

30

12

180

50

2

198

4

80

3.2

196.8

5

100

4

196

6

300

12

188

7

500

20

180

Diluted
1
2
3

√
√

Enough of each standard was prepared to load 5μL in triplicate per microplate.

107

Appendix I: Information needed to submit RNA for library construction and sequencing
to University of Illinois Keck Center.
* a word or excel file with: sample ID, label on tube, concentration (Qubit) and volume. Also, information
on how the samples will be pooled and sequenced (i.e. how many lanes, which samples should be pooled
in each lane, single-read or paired-end, etc). This form can be filled out completely to satisfy this request.
* Submit at least 1ug of total RNA in a minimum volume of 20ul and a maximum volume of 50ul of RNAsefree water
* total RNA should be free of genomic DNA (either DNAsed or otherwise free of contamination with
gDNA).
* run an aliquot (~50 to 100ng) of the total RNA on a 1% agarose gel next to a DNA ladder and send us a
picture before you submit the sample. This picture will be used to evaluate integrity of the total RNA as
well as presence/absence of gDNA. See representative gels and bioanalyzer traces
here: http://biotech.illinois.edu/htdna/services-equipment/illumina
* alternatively, run the total RNA on an RNA bioanalyzer chip or AATI Fragment Analyzer and send us the
.pdf file.
* on-campus only: account number to be used for the project. Charges will not be applied until after the
data is delivered but the account is needed to set up the project.
* off-campus only: Off-campus user form or TTA must be in place. Charges will be applied at contract
execution or after work is completed, depending upon negotiated terms.
*Shipping Address:
University of Illinois Keck Center
Attn: Chris Wright
1201 W. Gregory Dr.
334 ERML
Urbana, IL 61801
(217) 333-4372

1)
2)
3)
4)

PI Name: John M. Hranitz
Customer Name: Bloomsburg University of Pennsylvania
Date: April 8, 2019
FOAPAL (on-campus customers only):

108

Sample ID
(this name
will be used
to name you
fastq files).
Use only
letters,
numbers and
underscores.
Neg. 0-13
Neg. 0-14
Neg. 0-8
Neg. 0-11
Neg. 0-4
Neg. 0-5
Neg. 4-6
Neg. 4-7
Neg. 4-9
Neg. 4-10
Neg. 4-11
Neg. 4-5
Treat. 0-2
Treat. 0-3
Treat. 0-4
Treat. 0-7
Treat. 0-9
Treat. 0-6
Treat. 4-6
Treat. 4-8
Treat. 4-1
Treat. 4-5
Treat. 4-9
Treat. 4-3

Label on Tube
(we strongly
suggest #1-xx)

Concentration
(Qubit or
Nanodrop?), ng/ul

Volume, ul

Total RNA
(ng)

Final Pool #

I-A
I-B
I-C
I-D
I-E
I-F
II-A
II-B
II-C
II-D
II-E
II-F
III-A
III-B
III-C
III-D
III-E
III-F
IV-A
IV-B
IV-C
IV-D
IV-E
IV-F

2
1.92
1.86
1.86
1.85
1.8
1.73
1.8
1.66
1.66
1.66
1.64
1.78
2.05
1.94
1.93
1.89
1.8
1.9
1.91
1.87
1.84
1.83
1.8

30
30
30
30
30
30
30
30
30
30
30
30
30
30
30
30
30
30
30
30
30
30
30
30

35.8
36.2
45.4
34.8
36.5
39.2
47
56.1
48.6
46.5
31.2
24.5
50.7
47.9
56.8
57.9
40.9
34.7
24.9
25.2
50.3
39.9
22
34.7

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

5) Quality check of RNA: run an aliquot (~50 to 100ng) of the total RNA on a 1%
agarose gel next to a DNA ladder and attach or embed the image here. This picture
will be used to evaluate integrity of the total RNA as well as presence/absence of
gDNA. Alternatively, run the total RNA on an RNA bioanalyzer chip and attach the
.pdf file.
6) How many lanes for sequencing, type of run (single-read or paired-end) and length of
run (100nt, 2x150nt, etc):
*Number of lanes: 2
*Single or Paired-reads: Single
109

*Run type and read length: HiSeq 4000 100nt single-end read
MiSeq (1 lane per run):
Nano 2x250nt (500k – 1M paired reads)
Bulk 2x250nt (10-20M paired reads)
Bulk 2x300nt (25-50M paired reads)
NovaSeq options:

Flowcell

Read Type and
Length
Singlereads, 100bp
Paired-reads,
2x150bp

Price per # of Single- # of PairedLane
Reads
Reads
Gb/lane
$1,720

380
million

40

$2,720

750
million

120

Paired-reads,
2x250bp

$3,770

750
million

200

Singlereads, 100bp

$3,400

Pairedreads, 2x100bp

$4,890

1.5 billion

160

Paired-reads,
2x150bp

$5,270

1.5 billion

240

SP

750
million

80

S1

S2

Singlereads, 100bp

$6,450

S4

Paired-reads,
2x150bp

$8,990

1.5
billion

7) Any additional comments, special instructions, etc:

110

165

5-6 billion

750850

Appendix J: RNA-Seq Report prepared by the University of Illinois.
John M Hranitz RNA-Seq Report
Negin Valizadegan and Jenny Drnevich, HPCBio, University of Illinois
Jun 7, 2019
Contents
Location of results and codes for reproducibility .......................................................................................... 1
Quality Check, alignment and count generation ............................................................................................ 1
Basic statistical analysis ................................................................................................................................ 2
Hardware and software descriptions ............................................................................................................. 5
R and package versions ................................................................................................................................. 5
Location of results and codes for reproducibility
All deliverables for the basic analysis results are in a zipped file that can be downloaded from
box.com. Please unzip the file after transferring to your computer (on most PCs: right click and
select
“Extract
all”).
The
individual
gene
results
are
in
the
files
“jmhranitz_results_24Samples_2019-05-30.xlsx”, the raw gene-level counts (summed from all
transcripts) are in the file “2019-05-17_Gene_level_counts.xlsx” and the normalized logCPM values
for all samples are in the file “jmhranitz_gene_logCPMvalues_2019-05-30.xlsx”. Note that the .html
files in the “interactive_results” folder need to be kept with the “css” and “js” subfolders. The
“final_report” folder contains the file “Reports_jmhranitz_Honeybee_2019Apr.Rmd” that generated
this report and includes the R codes for the entire analysis. All necessary files to run the analysis,
including the raw transcript-level counts in .RData format are also in the “final_report” folder. The
.Rmd file was rendered in .html and pdf formats and the figures generated in each are also available
in the “Reports_jmhranitz_Honeybee_2019Apr_files” folder.
Quality Check, alignment and count generation
Reference sequences
The Apis mellifera transcriptome and Annotation Release 104 from NCBI are used for quasi-mapping
and count generation. This transcriptome is derived from genome Amel_HAv3.1. Since the quasimapping step only uses transcript sequences, the annotation file was solely used to generate transcriptgene mapping table which was kept in RData file “txID2GeneID.RData” for obtaining gene-level
counts.
Quality check on the raw data
We reviewed the QC report, “Project_rmhranitz_24_RNAseq_multiqc_report.html” sent by the
sequencing center that performed FASTQC 1 (version 0.11.8) on individual samples then was
summarized into a single html report by using MultiQC 2 version 1.6. Average per-base read quality
scores are over 30 in all samples and no adapter sequences were found indicating those reads are
high in quality. Thus, we skipped the trimming step and directly proceed to transcripts mapping
and quantification.
1
Andrews, S. (2010). FastQC: a quality control tool for high throughput sequence data. Available at:
https://www. bioinformatics.babraham.ac.uk/projects/fastqc/
2

Ewels, P., Magnusson, M., Lundin, S., & K?ller, M. (2016). MultiQC: summarize analysis results for
multiple tools and samples in a single report. Bioinformatics, 32(19), 3047-3048.

111

Percentage of reads mapped to transcriptome
100

Sam ple
90

80

70

Percentage

60

50
40

30

20

10

Figure 1: Figure 1. ReadFate plot
Alignment and gene-level quantification
Salmon3 version 0.13.1 was used to quasi-map reads to the transcriptome and quantify the abundance
of each transcript. The transcriptome was first indexed, then quasi-mapping was performed to map
reads to transcriptome with additional arguments --seqBias and --gcBias to correct sequence-specific
and GC content biases and --numBootstraps=30 to compute bootstrap transcript abundance
estimates. Gene-level counts were then estimated based on transcript-level counts using the “bias
corrected counts without an offset” method from the tximport package. This method provides more
accurate gene-level counts estimates and keeps multi-mapped reads in the analysis compared to
traditional alignment-based method4.
Percentage of reads mapped to the transcriptome ranged from 70.1 to 82.3% (Figure 1). The unmapped
reads were discarded while the number of remaining reads (range: 28.7 - 50.8 million per sample) were
kept for statistical analysis.

Basic statistical
analysis
Normalization and
filtering
When comparing expression levels, the numbers of reads per gene need to be normalized not only
because of the differences in total number of reads, but because there could be differences in RNA
112

Treat_4_9

Treat_4_8

Treat_4_6

Treat_4_5

Treat_4_3

Treat_4_1

Treat_0_9

Treat_0_7

Treat_0_6

Treat_0_4

Treat_0_3

Treat_0_2

Neg_4_9

Neg_4_7

Neg_4_6

Neg_4_5

Neg_4_11

Neg_4_10

Neg_0_8

Neg_0_5

Neg_0_4

Neg_0_14

Neg_0_13

Neg_0_11

0

composition such that the total number of reads would not be expected to be the same. The TMM
(trimmed mean of M
3
Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., & Kingsford, C. (2017). Salmon provides fast and bias-aware quantification
of transcript expression. Nature methods, 14(4), 417.
4
Soneson C, Love MI, Robinson MD (2015). “Differential analyses for RNA-seq: transcript-level estimates improve genelevel inferences.” F1000Research, 4. doi: 10.12688/f1000research.7563.1.

113

values) normalization5 in the edgeR package6 uses the assumption of most genes do not change to calculate
a normalization factor for each sample to adjust for such biases in RNA composition. In this dataset,
TMM normalization factors fluctuates between 0.81 and 1.2 but the variation is between individuals,
not groups, which suggests no group-level overall differences in RNA composition. While the NCBI
Amel_HAv3.1 Annotation Release 104 gene models have a total of 12,090 genes, some of these might
not have detectable expression. Therefore, we set the detection threshold at 0.5 cpm (counts per
million) in at least 3 samples, which resulted in 1,493 genes being filtered out, leaving 10,597 genes to
be analyzed for differential expression that contain 99.99% of the reads. After filtering, TMM
normalization was performed again and normalized log2-based count per million values (logCPM)
were calculated using edgeR's cpm() function with prior.count = 3 to help stabilize fold-changes of
extremely low expression genes.
Clustering
Multidimensional scaling in the limma7 package was used to identify potential treatment effects at higher
level. The normalized logCPM values of the top 5,000 variable genes were chosen to construct the
multidimensional scaling plot (Figure 3). Dimension 1 which explains around —25% of the total
variability, explains the separation of negative control at time 4 from the rest of the groups. Dimension
2 which contains —17% of the variability does not explain groupings of the effect of treatment, negative
control and/or time. An interactive version of this plot can be found at
“interactive_results/MDSclustering_postFiltering.html”.
Robinson MD, Oshlack A (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome
Biology 11, R25.

5

Robinson MD, McCarthy DJ, Smyth GK (2010). “edgeR: a Bioconductor package for differential expression analysis of digital gene
expression data.” Bioinformatics, 26(1), 139-140.

6

114

7
Ritchie, M.E., Phipson, B., Wu, D., Hu, Y., Law, C.W., Shi, W., and Smyth, G.K. (2015). limma powers differential
expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 43(7), e47

115

116

Table 1: Table 1. Number of differenially expressed genes (global FDR p < 0.1)
Down
NotSig
Up

T4vsT0

N4vsN0

T0vsN0

T4vsN4

Interact

46

2369

35

1663

1091

10476

6312

10365

6725

7805

75

1916

197

2209

1701

Differential expression testing
Differential gene expression (DE) analysis was performed using the limma-trend method8 9 and all
four logical pairwise comparisons were computed, along with a test for any interaction between
treatment and time. Because pairwise comparisons involving the Neg_4 group had so many more
DE genes than comparisons between other groups, we adjust for multiple testing correction by
doing a “global” False Discovery Rate correction10 across p-values for all 5 comparisons together.
This ensures that a gene with the same raw p-value in two different comparisons would not end
up with vastly different FDR p-values. The number of up and down regulated genes using FDR pvalue < 0.1 for the four logical pairwise comparisons and the interaction test are in Table 1.
Interactive versions of the results for each comparison are in “.html” files in the
“interactive_results” folder. We also calculated the equivalent of a oneway ANOVA test across all
four groups to select genes to visualize their overall expression patterns in a heatmap. False
discovery rate correction was done separately for the oneway ANOVA, and at FDR < 0.05, 3819
genes showed differences across the groups, mainly the Neg_4 group was different from the others
(Figure 4).
Functional annotation
Gene name, Gene Ontology ID and GO terms for each gene was obtained by using the
org.Apis_mellifera.eg.sqlite database from AnnotationHub package from Bioconductor11 release 3.9.
KEGG pathways for each gene were retrieved directly from http://www.genome.jp/kegg using the
KEGGREST package on May 30, 2019 .

Hardware and software descriptions
The read quality check and count generation was done using CNRG's Biocluster high-performance
computing resource. All analyses from summation of counts to the gene level (including the codes
in .Rmd file) were done on a laptop in R version 3.6.0 (2019-04-26)12 using packages as indicated
below.

R and package versions
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
8
Chen Y, Lun ATL and Smyth GK. From reads to genes to pathways: differential expression analysis of RNA-Seq
experiments using Rsubread and the edgeR quasi-likelihood pipeline [version 2; referees: 5 approved].
F1000Research 2016, 5:1438 (doi: 10.12688/f1000research.8987.2)
9
Law, C. W., Chen, Y., Shi, W., & Smyth, G. K. (2014). voom: Precision weights unlock linear model analysis tools
for RNA-seq read counts. Genome biology, 15(2), R29.
10
Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to
multiple testing. Journal of the Royal statistical society: series B (Methodological), 57(1), 289-300.

117

Gentleman, R. C., Carey, V. J., Bates, D. M., Bolstad, B., Dettling, M., Dudoit, S., ... & Hornik, K. (2004).
Bioconductor: open software development for computational biology and bioinformatics. Genome biology, 5(10),
R80.

11

12R Core Team (2017). R: A language and environment for statistical computing. R Foundation for Statistical

Computing, Vienna, Austria. URL https://www.R-project.org/.

118

−2 0 2

119

Neg_4_5
Neg_4_6
Neg_4_7
Neg_4_9
Treat_0_2
Treat_0_3
Treat_0_4
Treat_0_6
Treat_0_7
Treat_0_9
Treat_4_1
Treat_4_3
Treat_4_5
Treat_4_6
Treat_4_8

Neg_0_11
Neg_0_13
Neg_0_14
Neg_0_4
Neg_0_5
Neg_0_8
Neg_4_10
Neg_4_11

Color Key

3819 genes with
1way ANOVA FDR p < 0.05

SD from mean

Figure 4: Figure 4. Heatmap of expression patterns of genes with a significant one-way ANOVA

Appendix K: Sequence quality from Multi-QC report.

120

Appendix L: Outlier genes removed from each comparison group.
Control 4 HourTreatment 0 Hour

Treatment 4
Hour-Control 0
Hour

Treatment 0
Hour-Control 0
Hour

Treatment 4
Hour-Control 4
Hour

LOC102654530
LOC102655462
LOC408310
LOC726513
LOC100578618
LOC102656697
Or63
LOC102654769
LOC725762
LOC100576417
LOC100577069
LOC102656780
LOC102656729
LOC102655825
LOC100577908

LOC100577041
LOC107964312
LOC113219186
LOC113218895
Mir6058
LOC107965500
LOC102654676
LOC100578137
LOC102655174
LOC102655148
LOC113218910
LOC100578882
LOC100576914
LOC413532
LOC113219427
LOC102654659
LOC413574
LOC727035
LOC100577278
LOC102653680
LOC102656189
Mir3726
LOC102655107
LOC102653600
LOC102654231
LOC107965601
LOC102656705
LOC102653975
LOC107965171
LOC100577817
LOC100578397
Mir9865
LOC725965
LOC102654917
LOC100577474
LOC100577042
LOC113219125
LOC725469
Or25
LOC102655232
Mir6049
LOC726468
LOC726850
LOC107964633
LOC100576944
Mir125
LOC107965228
LOC113219070
LOC102654757
LOC113219128
LOC102656145
LOC113219009
LOC725606
Eyg
LOC102654300

LOC102654769
LOC102656780
LOC107965596
LOC107965494
LOC100577266
LOC102654837
LOC100577278
LOC113218598
LOC100576599
LOC113219192
LOC100577149
LOC102654300
LOC107965117
Or63
LOC102656697
LOC102654810
LOC102656919
LOC107964750
LOC102656729
LOC102654804
LOC107964633
Mir9883
LOC102654490
Apd-2
LOC410603

LOC113218730
LOC727126
LOC107964028
LOC113218895
LOC102654676
LOC725965
LOC100577041
LOC107966041
LOC107965500
LOC725469
LOC410107
LOC413166
LOC102655148
LOC107964312
LOC725762
LOC413574
Mir6042
LOC113219059
LOC107964934
LOC727035
LOC102655107
LOC102655174
LOC107964975
LOC724312
Mir6058
LOC100576914
LOC113218527
LOC727335
LOC102654660
LOC724895
LOC100578446
LOC113218978
LOC727642
LOC102654659
LOC100577521
LOC100578137
LOC113218765
LOC113219186
LOC551263
LOC726072
LOC726863
LOC113219196
Or105
LOC102656294
LOC107964152
LOC724860
LOC409751
LOC410231
LOC724386
LOC100577752
LOC100576082
LOC102655724
LOC102656479
LOC102656705
LOC102653680

121

Mir3759
LOC102654929
LOC102655235
LOC100577149
LOC102654924
LOC724895
LOC113219151
LOC727642
LOC113218881
LOC113219388
LOC107964043
LOC113218774
LOC410687
LOC102653852
LOC107965367
LOC100578057
LOC100577559
LOC107964737
LOC107965375
LOC113218672
LOC107965480
LOC107964483
LOC100577266
LOC724429
LOC100578177

LOC107964737
LOC113219200
LOC107965534
LOC113219050
LOC102653963
LOC113218668
LOC113218910
LOC725945
LOC100576688
LOC107965073
LOC100576172
LOC410254
LOC100578882
Mir6049
LOC102656719
LOC100578177
LOC724645
Mir3726
LOC107964006
LOC102653852
LOC100578193
LOC413532
LOC107965228
LOC113219340
LOC724934
Mir9865
LOC113218946
LOC113219185
LOC725831
LOC113218795
LOC724921
LOC113218522
LOC113219051
LOC113219388
LOC113219370
LOC725935
LOC100576901
LOC100576609
LOC724624
LOC113218653
LOC102654144
LOC725589
LOC102655670
LOC102656570
LOC100576944
LOC102655155
LOC107965158
LOC102656509
LOC100576439
LOC107965805
LOC102653600
LOC107965171
LOC410334
LOC100578130
LOC100577069
LOC113219009
LOC113219128
Eyg
Mir125
LOC413693
LOC107965375
LOC113218881

122

LOC410736
LOC727237
LOC100577817
LOC724429
LOC102654622
LOC102654948
LOC725302
LOC102656864
LOC411387
LOC100577373
LOC100578474
LOC410687
LOC113218672
LOC113219125
LOC113219070
LOC100577559
LOC102656387
LOC100576387
LOC551632
LOC100578729

123

Appendix M: Python script.

124

Appendix N: Average Motor Scores of Honey Bees Collected and Treated with
Imidacloprid

Treatment
Groups

Concentration
of Imidacloprid
(ng/10 µL)

Average
Abdomen
Score

Average
leg score

Average
Antennae
Score

Average
P.E.R.
score

Average
Motor
Score

Number
of Bees
Treated

Negative

0

1.75

1.95

2.00

2.00

7.70

20

1/5

3.6

1.26

1.74

1.53

1.58

6.11

19

1/10

1.8

1.67

1.78

1.78

1.61

6.83

18

1/20

0.9

1.82

1.88

1.82

1.88

7.41

17

1/50

0.36

1.70

1.90

1.95

1.90

7.45

20

1/100

0.18

1.65

1.95

2.00

1.95

7.55

20

1/500

0.036

1.82

1.71

2.00

1.94

7.47

17

Positive

0

1.85

1.69

1.69

1.77

7.00

13

125

Appendix O: List of 100 most significant changes in Treatment-4h versus Control-0h.

126

kelch domain-containing protein 10 homolog

guanine nucleotide exchange factor for Rab-3A, transcript variant X1

vesicle-fusing ATPase 1

ethanolamine-phosphate cytidylyltransferase, transcript variant X1

serine/threonine-protein kinase 3, transcript variant X1

staphylococcal nuclease domain-containing protein 1

phosphatidylinositol 5-phosphate 4-kinase type-2 alpha, transcript variant X2

ribonuclease Z, mitochondrial, transcript variant X1

DCN1-like protein 3

E3 ubiquitin-protein ligase ZNRF2

fizzy-related protein homolog

sulfhydryl oxidase 1

WD repeat and FYVE domain-containing protein 2

potential E3 ubiquitin-protein ligase ariadne-1, transcript variant X3

dnaJ homolog subfamily C member 16

dnaJ homolog subfamily C member 3

homer protein homolog 2, transcript variant X1

uncharacterized LOC413653, transcript variant X2

ubiquitin carboxyl-terminal hydrolase 14

uncharacterized LOC100576438, transcript variant X2

sorting nexin lst-4

myb-like protein D

CTD small phosphatase-like protein 2, transcript variant X3

methyltransferase-like protein 23

uncharacterized LOC412397, transcript variant X2

calmodulin

rho GTPase-activating protein 1, transcript variant X1

arfaptin-2, transcript variant X2

serine/threonine-protein phosphatase 2A catalytic subunit alpha isoform

cyclin-Y
RING finger and CHY zinc finger domain-containing protein 1,
transcript variant X3

tubulin beta chain
lysophospholipase D GDPD1-like, transcript variant X4

galactoside 2-alpha-L-fucosyltransferase 2

dual specificity mitogen-activated protein kinase kinase dSOR1

beta-1,3-galactosyltransferase 1-like, transcript variant X2

syntaxin-12, transcript variant X3
actin-related protein 2, transcript variant X1

transmembrane protein 64
acidic leucine-rich nuclear phosphoprotein 32 family member A,
transcript variant X2

peroxisomal N(1)-acetyl-spermine/spermidine oxidase

pantothenate kinase 3, transcript variant X1

death-associated inhibitor of apoptosis 2, transcript variant X2

WD repeat-containing protein 82, transcript variant X1

leucine-rich repeat protein soc-2 homolog

protein HID1, transcript variant X2

uncharacterized LOC410622

WD repeat-containing protein 48

transmembrane protein 62, transcript variant X1

DNA fragmentation factor subunit alpha

126

MAPK regulated corepressor interacting protein 2

hyccin, transcript variant X3

BTB/POZ domain-containing adapter for CUL3-mediated RhoA degradation protein CUE
3, transcript
domain-containing
variant X2 protein 2

127

N-alpha-acetyltransferase 30

blastoderm-specific protein 25D

transmembrane 9 superfamily member 3

zinc finger BED domain-containing protein 1

UV excision repair protein RAD23 homolog B, transcript variant X2

ras-related protein Rab-30

ras association domain-containing protein 2

NECAP-like protein CG9132, transcript variant X1

GTP-binding protein Rit2, transcript variant X3

WD repeat domain phosphoinositide-interacting protein 3

alpha-1,3/1,6-mannosyltransferase ALG2

uncharacterized LOC552002, transcript variant X3

zwei Ig domain protein zig-8, transcript variant X2

RILP-like protein homolog, transcript variant X2

choline-phosphate cytidylyltransferase A, transcript variant X3

probable galactose-1-phosphate uridylyltransferase

dnaJ homolog subfamily C member 7, transcript variant X2

probable E3 ubiquitin-protein ligase makorin-1, transcript variant X1

protein zyg-11 homolog B, transcript variant X4

ethanolaminephosphotransferase 1

V-type proton ATPase subunit C, transcript variant X4

cleft lip and palate transmembrane protein 1 homolog

uncharacterized LOC724843, transcript variant X2

NADPH--cytochrome P450 reductase, transcript variant X2

uncharacterized LOC551031, transcript variant X1

uncharacterized LOC408732

ceramide synthase 6, transcript variant X1

condensin complex subunit 2

ER degradation-enhancing alpha-mannosidase-like protein 2

adipocyte plasma membrane-associated protein

protein spinster, transcript variant X2

calreticulin

actin-binding protein IPP

AFG3-like protein 2

serine/threonine-protein kinase PINK1, mitochondrial, transcript variant X2

E3 ubiquitin-protein ligase NRDP1, transcript variant X2

DCN1-like protein 5, transcript variant X1

protein disulfide-isomerase

palmitoyltransferase app, transcript variant X2

SH3 domain-binding glutamic acid-rich protein homolog, transcript variant X2

dual specificity protein phosphatase CDC14AB, transcript variant X3

O-phosphoseryl-tRNA(Sec) selenium transferase, transcript variant X3

mitoferrin-1, transcript variant X1

zinc finger protein 25

TOM1-like protein 2, transcript variant X1

WD repeat domain phosphoinositide-interacting protein 2, transcript variant X2

127

Appendix P: Control Time 4 Hour-Treatment Time 0 Hour KEGG pathways. Genes that are colored red are considered downregulated (not significant), genes that are colored purple are significantly down-regulated, genes that are colored pink are significantly
up-regulated, and genes that are colored orange are considered up-regulated (not significant) (significance p<0.1). Enzymes are
indicated by the green color.

128
128

129
129

130
130

131
131

132
132

133
133

134
134

135
135

136
136

137
137

Appendix Q: Treatment Time 4 Hour-Control Time 0 Hour KEGG pathways. Genes that are colored red are considered downregulated (not significant), genes that are colored purple are significantly down-regulated, genes that are colored pink are significantly
up-regulated, and genes that are colored orange are considered up-regulated (not significant) (significance p<0.1). Enzymes are
indicated by the green color.

138
138

139
139

140
140

141
141

142
142

143
143

144
144

145
145

Appendix R: DAVID Analysis of significantly (p<0.1) enriched genes for Treatment-4h versus
Control-0h.
Down-Regulated
Annotation Cluster 1
Category
GOTERM_MF_DIRECT
INTERPRO
INTERPRO
SMART

Enrichment Score: 5.45637131589216
Term
GO:0008270~zinc ion binding
IPR013083:Zinc finger, RING/FYVE/PHD-type
IPR001841:Zinc finger, RING-type
SM00184:RING

Annotation Cluster 2
Category
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS
GOTERM_CC_DIRECT

Enrichment Score: 4.8010075530818375
Term
Membrane
Transmembrane helix
Transmembrane
GO:0016021~integral component of membrane

Annotation Cluster 3
Category
GOTERM_BP_DIRECT
INTERPRO
INTERPRO
GOTERM_MF_DIRECT
UP_KEYWORDS
INTERPRO

Enrichment Score: 2.942701935144724
Term
GO:0007264~small GTPase mediated signal transduction
IPR005225:Small GTP-binding protein domain
IPR001806:Small GTPase superfamily
GO:0005525~GTP binding
GTP-binding
IPR027417:P-loop containing nucleoside triphosphate hydrolase

Annotation Cluster 4
Category
GOTERM_MF_DIRECT
UP_KEYWORDS
GOTERM_MF_DIRECT
INTERPRO
SMART
SMART

Enrichment Score: 2.800889694997264
Term
GO:0004842~ubiquitin-protein transferase activity
Ubl conjugation pathway
GO:0016874~ligase activity
IPR000569:HECT
SM00119:HECTc
SM00456:WW

Annotation Cluster 5
Category
INTERPRO
INTERPRO
UP_KEYWORDS
UP_KEYWORDS
GOTERM_MF_DIRECT
INTERPRO
SMART
INTERPRO
UP_KEYWORDS
UP_KEYWORDS
GOTERM_MF_DIRECT

Enrichment Score: 2.7931918200278925
Term
IPR008271:Serine/threonine-protein kinase, active site
IPR017441:Protein kinase, ATP binding site
Kinase
Serine/threonine-protein kinase
GO:0004674~protein serine/threonine kinase activity
IPR000719:Protein kinase, catalytic domain
SM00220:S_TKc
IPR011009:Protein kinase-like domain
Nucleotide-binding
ATP-binding
GO:0005524~ATP binding

146

Annotation Cluster 6
Category
GOTERM_BP_DIRECT
INTERPRO
SMART

Enrichment Score: 2.777276800524685
Term
GO:0035556~intracellular signal transduction

Annotation Cluster 7
Category
INTERPRO
INTERPRO
SMART

Enrichment Score: 1.9311441030506946
Term
IPR011993:Pleckstrin homology-like domain
IPR001849:Pleckstrin homology domain
SM00233:PH

Annotation Cluster 8
Category
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 1.7901506620531265
Term
Zinc
Metal-binding
Zinc-finger

Annotation Cluster 9
Category
INTERPRO
INTERPRO
SMART

Enrichment Score: 1.579437655647095
Term
IPR011333:BTB/POZ fold
IPR000210:BTB/POZ-like
SM00225:BTB

Annotation Cluster 10
Category
SMART
SMART
INTERPRO

Enrichment Score: 1.571201116530967
Term
SM00253:SOCS
SM00969:SM00969
IPR001496:SOCS protein, C-terminal

Annotation Cluster 11
Category
INTERPRO
SMART
UP_KEYWORDS

Enrichment Score: 1.4532259278548885
Term
IPR001452:Src homology-3 domain
SM00326:SH3
SH3 domain

Annotation Cluster 12
Category
INTERPRO
INTERPRO
INTERPRO
INTERPRO
INTERPRO
GOTERM_MF_DIRECT
INTERPRO

Enrichment Score: 1.2000790546574989
Term
IPR023214:HAD-like domain
IPR023299:P-type ATPase, cytoplasmic domain N
IPR001757:Cation-transporting P-type ATPase
IPR018303:P-type ATPase, phosphorylation site
IPR008250:P-type ATPase, A domain
GO:0004012~phospholipid-translocating ATPase activity
IPR006539:Phospholipid-transporting P-type ATPase, subfamily IV

Annotation Cluster 13
Category

Enrichment Score: 1.191896405979736
Term

SM00109:C1

147

INTERPRO
GOTERM_MF_DIRECT
INTERPRO
INTERPRO
SMART

IPR011992:EF-hand-like domain
GO:0005509~calcium ion binding
IPR002048:EF-hand domain
IPR018247:EF-Hand 1, calcium-binding site
SM00054:EFh

Annotation Cluster 14
Category
INTERPRO
INTERPRO
UP_KEYWORDS
SMART
INTERPRO
SMART

Enrichment Score: 1.1422749230721754
Term
IPR013763:Cyclin-like
IPR006671:Cyclin, N-terminal
Cyclin
SM00385:CYCLIN
IPR004367:Cyclin, C-terminal domain
SM01332:SM01332

Annotation Cluster 15
Category
INTERPRO
INTERPRO
SMART

Enrichment Score: 1.0468816164360077
Term
IPR017892:Protein kinase, C-terminal
IPR000961:AGC-kinase, C-terminal
SM00133:S_TK_X

Annotation Cluster 16
Category
INTERPRO
SMART
INTERPRO

Enrichment Score: 1.044975678587419
Term
IPR013761:Sterile alpha motif/pointed domain
SM00454:SAM
IPR001660:Sterile alpha motif domain

Annotation Cluster 17
Category
UP_KEYWORDS
INTERPRO
SMART
SMART
INTERPRO
INTERPRO
INTERPRO
GOTERM_MF_DIRECT

Enrichment Score: 0.9691442961294987
Term
Protein phosphatase
IPR003595:Protein-tyrosine phosphatase, catalytic
SM00404:PTPc_motif
SM00194:PTPc

Annotation Cluster 18
Category
INTERPRO
GOTERM_BP_DIRECT
INTERPRO
GOTERM_BP_DIRECT
GOTERM_MF_DIRECT
INTERPRO

Enrichment Score: 0.9629833272596152
Term
IPR001394:Peptidase C19, ubiquitin carboxyl-terminal hydrolase 2
GO:0016579~protein deubiquitination

Annotation Cluster 19
Category

Enrichment Score: 0.9218108166394164
Term

IPR016130:Protein-tyrosine phosphatase, active site
IPR000387:Protein-tyrosine/Dual specificity phosphatase
GO:0004725~protein tyrosine phosphatase activity

GO:0006511~ubiquitin-dependent protein catabolic process
GO:0036459~thiol-dependent ubiquitinyl hydrolase activity
IPR001607:Zinc finger, UBP-type

148

INTERPRO
SMART
INTERPRO
GOTERM_BP_DIRECT

IPR000300:Inositol polyphosphate-related phosphatase
SM00128:IPPc
IPR005135:Endonuclease/exonuclease/phosphatase
GO:0046856~phosphatidylinositol dephosphorylation

Annotation Cluster 20
Category
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.9059256645404106
Term
IPR018108:Mitochondrial substrate/solute carrier
IPR023395:Mitochondrial carrier domain
IPR002067:Mitochondrial carrier protein

Annotation Cluster 21
Category
UP_KEYWORDS
INTERPRO
SMART

Enrichment Score: 0.8667189123896825
Term
LIM domain
IPR001781:Zinc finger, LIM-type
SM00132:LIM

Annotation Cluster 22
Category
KEGG_PATHWAY
KEGG_PATHWAY
KEGG_PATHWAY

Enrichment Score: 0.8600660685290529
Term
ame00072:Synthesis and degradation of ketone bodies
ame00650:Butanoate metabolism
ame00280:Valine, leucine and isoleucine degradation

Annotation Cluster 23
Category
GOTERM_BP_DIRECT
INTERPRO
GOTERM_CC_DIRECT

Enrichment Score: 0.8233643535747749
Term
GO:0006888~ER to Golgi vesicle-mediated transport
IPR006895:Zinc finger, Sec23/Sec24-type
GO:0030127~COPII vesicle coat

Annotation Cluster 24
Category
GOTERM_BP_DIRECT
GOTERM_CC_DIRECT
GOTERM_MF_DIRECT

Enrichment Score: 0.8135565670707758
Term
GO:0015991~ATP hydrolysis coupled proton transport
GO:0033179~proton-transporting V-type ATPase, V0 domain
GO:0015078~hydrogen ion transmembrane transporter activity

Annotation Cluster 25
Category
KEGG_PATHWAY
GOTERM_MF_DIRECT
INTERPRO
INTERPRO
INTERPRO
SMART
GOTERM_MF_DIRECT
GOTERM_BP_DIRECT

Enrichment Score: 0.7964987710011121
Term
ame01210:2-Oxocarboxylic acid metabolism
GO:0004449~isocitrate dehydrogenase (NAD+) activity
IPR024084:Isopropylmalate dehydrogenase-like domain
IPR001804:Isocitrate and isopropylmalate dehydrogenases family
IPR004434:Isocitrate dehydrogenase NAD-dependent
SM01329:SM01329
GO:0051287~NAD binding
GO:0006099~tricarboxylic acid cycle

Annotation Cluster 26
Category
INTERPRO

Enrichment Score: 0.7867348565935095
Term
IPR001623:DnaJ domain

149

SMART
INTERPRO

SM00271:DnaJ
IPR018253:DnaJ domain, conserved site

Annotation Cluster 27
Category
INTERPRO
GOTERM_MF_DIRECT
INTERPRO
SMART

Enrichment Score: 0.7709951919258211
Term
IPR004843:Metallophosphoesterase domain
GO:0004721~phosphoprotein phosphatase activity

Annotation Cluster 28
Category
GOTERM_BP_DIRECT
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.7671408489363004
Term
GO:0015991~ATP hydrolysis coupled proton transport
IPR004100:ATPase, alpha/beta subunit, N-terminal

Annotation Cluster 29
Category
INTERPRO
UP_KEYWORDS
INTERPRO
INTERPRO
INTERPRO
GOTERM_BP_DIRECT
GOTERM_CC_DIRECT

Enrichment Score: 0.7477835262115027
Term
IPR017937:Thioredoxin, conserved site
Redox-active center
IPR005788:Disulphide isomerase
IPR013766:Thioredoxin domain
IPR012336:Thioredoxin-like fold
GO:0045454~cell redox homeostasis
GO:0005623~cell

Annotation Cluster 30
Category
INTERPRO
GOTERM_MF_DIRECT
SMART

Enrichment Score: 0.7103214249214863
Term
IPR001683:Phox homologous domain
GO:0035091~phosphatidylinositol binding
SM00312:PX

Annotation Cluster 31
Category
GOTERM_MF_DIRECT
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.6979412102539907
Term
GO:0016791~phosphatase activity
IPR011948:Dullard phosphatase domain, eukaryotic
IPR004274:NLI interacting factor
SM00577:CPDc

Annotation Cluster 32
Category
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.6930615590326878
Term
IPR010504:Arfaptin homology (AH) domain
SM01015:SM01015
IPR027267:Arfaptin homology (AH) domain/BAR domain

Annotation Cluster 33
Category
INTERPRO

Enrichment Score: 0.6577784257414996
Term
IPR000301:Tetraspanin

SM00156:PP2Ac

150

INTERPRO
PIR_SUPERFAMILY
INTERPRO
INTERPRO

IPR018499:Tetraspanin/Peripherin
PIRSF002419:tetraspanin
IPR018503:Tetraspanin, conserved site
IPR008952:Tetraspanin, EC2 domain

Annotation Cluster 34
Category
KEGG_PATHWAY
KEGG_PATHWAY
GOTERM_BP_DIRECT

Enrichment Score: 0.645183527714244
Term
ame01210:2-Oxocarboxylic acid metabolism
ame00220:Arginine biosynthesis
GO:0006520~cellular amino acid metabolic process

Annotation Cluster 35
Category
GOTERM_MF_DIRECT
SMART
INTERPRO
KEGG_PATHWAY

Enrichment Score: 0.6186519111303531
Term
GO:0016746~transferase activity, transferring acyl groups
SM00563:PlsC
IPR002123:Phospholipid/glycerol acyltransferase
ame00561:Glycerolipid metabolism

Annotation Cluster 36
Category
INTERPRO
SMART
GOTERM_MF_DIRECT
GOTERM_MF_DIRECT

Enrichment Score: 0.5641078769792714
Term
IPR004827:Basic-leucine zipper domain
SM00338:BRLZ

Annotation Cluster 37
Category
GOTERM_BP_DIRECT
UP_KEYWORDS
UP_KEYWORDS
GOTERM_MF_DIRECT
INTERPRO
GOTERM_CC_DIRECT

Enrichment Score: 0.5448123900070027
Term
GO:0006486~protein glycosylation
Golgi apparatus
Glycosyltransferase
GO:0008378~galactosyltransferase activity
IPR002659:Glycosyl transferase, family 31
GO:0000139~Golgi membrane

Annotation Cluster 38
Category
INTERPRO
INTERPRO
INTERPRO
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.5379325126117979
Term
IPR017984:Chromo domain subgroup
IPR023779:Chromo domain, conserved site
IPR000953:Chromo domain/shadow
IPR016197:Chromo domain-like
SM00298:CHROMO
IPR023780:Chromo domain

Annotation Cluster 39
Category
GOTERM_MF_DIRECT
INTERPRO
INTERPRO
GOTERM_MF_DIRECT

Enrichment Score: 0.5227346997205058
Term
GO:0004871~signal transducer activity
IPR011025:G protein alpha subunit, helical insertion

GO:0043565~sequence-specific DNA binding

GO:0003924~GTPase activity

151

Annotation Cluster 40
Category
SMART
INTERPRO
INTERPRO
KEGG_PATHWAY
GOTERM_BP_DIRECT

Enrichment Score: 0.5110314775368822
Term
SM00397:t_SNARE
IPR000727:Target SNARE coiled-coil domain
IPR010989:t-SNARE
ame04130:SNARE interactions in vesicular transport
GO:0061025~membrane fusion

Annotation Cluster 41
Category
INTERPRO
INTERPRO
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.4881309535421234
Term
IPR011011:Zinc finger, FYVE/PHD-type
IPR001965:Zinc finger, PHD-type
IPR019787:Zinc finger, PHD-finger
SM00249:PHD
IPR019786:Zinc finger, PHD-type, conserved site

Annotation Cluster 42
Category
INTERPRO
INTERPRO
GOTERM_MF_DIRECT

Enrichment Score: 0.4842903939578597
Term
IPR001375:Peptidase S9, prolyl oligopeptidase, catalytic domain
IPR002469:Peptidase S9B, dipeptidylpeptidase IV N-terminal
GO:0008236~serine-type peptidase activity

Annotation Cluster 43
Category
SMART
INTERPRO
INTERPRO

Enrichment Score: 0.4798650796978842
Term
SM00558:JmjC
IPR003347:JmjC domain
IPR019786:Zinc finger, PHD-type, conserved site

Annotation Cluster 44
Category
GOTERM_MF_DIRECT
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS
KEGG_PATHWAY

Enrichment Score: 0.46779828181778016
Term
GO:0000049~tRNA binding
Protein biosynthesis
tRNA-binding
Aminoacyl-tRNA synthetase
Ligase
ame00970:Aminoacyl-tRNA biosynthesis

Annotation Cluster 45
Category
INTERPRO
INTERPRO
GOTERM_BP_DIRECT
INTERPRO

Enrichment Score: 0.44565558536911126
Term
IPR017983:GPCR, family 2, secretin-like, conserved site
IPR000832:GPCR, family 2, secretin-like
GO:0007166~cell surface receptor signaling pathway
IPR017981:GPCR, family 2-like

Annotation Cluster 46
Category
UP_KEYWORDS
GOTERM_BP_DIRECT

Enrichment Score: 0.42963092505039946
Term
Autophagy
GO:0006914~autophagy

152

KEGG_PATHWAY

ame04140:Regulation of autophagy

Annotation Cluster 47
Category
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.42745931692292394
Term
IPR001611:Leucine-rich repeat
IPR003591:Leucine-rich repeat, typical subtype
SM00369:LRR_TYP

Annotation Cluster 48
Category
INTERPRO
INTERPRO
INTERPRO
INTERPRO
PIR_SUPERFAMILY
SMART
INTERPRO
UP_KEYWORDS
SMART

Enrichment Score: 0.411021043275444
Term
IPR011043:Galactose oxidase/kelch, beta-propeller
IPR006652:Kelch repeat type 1
IPR011705:BTB/Kelch-associated
IPR017096:Kelch-like protein, gigaxonin
PIRSF037037:kelch-like protein, gigaxonin type
SM00875:SM00875
IPR015916:Galactose oxidase, beta-propeller
Kelch repeat
SM00612:Kelch

Annotation Cluster 49
Category
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.3926687745465885
Term
IPR020084:NUDIX hydrolase, conserved site
IPR000086:NUDIX hydrolase domain
IPR015797:NUDIX hydrolase domain-like

Annotation Cluster 50
Category
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.39040159347772047
Term
IPR017455:Zinc finger, FYVE-related
IPR000306:Zinc finger, FYVE-type
SM00064:FYVE

Annotation Cluster 51
Category
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.3892603002223318
Term
IPR013017:NHL repeat, subgroup
IPR001258:NHL repeat
IPR011042:Six-bladed beta-propeller, TolB-like

Annotation Cluster 52
Category
UP_KEYWORDS
GOTERM_MF_DIRECT
INTERPRO

Enrichment Score: 0.3805730839438125
Term
Acyltransferase
GO:0019706~protein-cysteine S-palmitoyltransferase activity
IPR001594:Zinc finger, DHHC-type, palmitoyltransferase

Annotation Cluster 53
Category
INTERPRO
UP_KEYWORDS
INTERPRO

Enrichment Score: 0.37907219453909835
Term
Pyridoxal phosphate
IPR015424:Pyridoxal phosphate-dependent transferase

153

INTERPRO
GOTERM_MF_DIRECT
INTERPRO

IPR004839:Aminotransferase, class I/classII
GO:0030170~pyridoxal phosphate binding

Annotation Cluster 54
Category
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS
GOTERM_MF_DIRECT
GOTERM_MF_DIRECT
GOTERM_MF_DIRECT
INTERPRO
UP_KEYWORDS
GOTERM_BP_DIRECT

Enrichment Score: 0.3747676070971319
Term
Lipid biosynthesis
Fatty acid metabolism
Fatty acid biosynthesis
GO:0102338~3-oxo-lignoceronyl-CoA synthase activity
GO:0102337~3-oxo-cerotoyl-CoA synthase activity
GO:0102336~3-oxo-arachidoyl-CoA synthase activity
IPR002076:GNS1/SUR4 membrane protein
Lipid metabolism
GO:0006633~fatty acid biosynthetic process

Annotation Cluster 55
Category
KEGG_PATHWAY
INTERPRO
GOTERM_BP_DIRECT
GOTERM_MF_DIRECT
KEGG_PATHWAY

Enrichment Score: 0.3710498061223294
Term
ame00220:Arginine biosynthesis
IPR004839:Aminotransferase, class I/classII
GO:0009058~biosynthetic process
GO:0030170~pyridoxal phosphate binding
ame00250:Alanine, aspartate and glutamate metabolism

Annotation Cluster 56
Category
GOTERM_MF_DIRECT
UP_KEYWORDS
INTERPRO
UP_KEYWORDS
GOTERM_CC_DIRECT
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 0.3592894436076364
Term
GO:0016805~dipeptidase activity
Dipeptidase
IPR008257:Peptidase M19, renal dipeptidase
GPI-anchor
GO:0031225~anchored component of membrane
Lipoprotein
Metalloprotease
Glycoprotein

Annotation Cluster 57
Category
SMART
INTERPRO
INTERPRO

Enrichment Score: 0.3556870936584284
Term
SM00060:FN3
IPR003961:Fibronectin, type III
IPR013783:Immunoglobulin-like fold

Annotation Cluster 58
Category
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.34945909361066296
Term
IPR004046:Glutathione S-transferase, C-terminal
IPR010987:Glutathione S-transferase, C-terminal-like
IPR004045:Glutathione S-transferase, N-terminal

Annotation Cluster 59
Category

Enrichment Score: 0.34163801579538106
Term

154

UP_KEYWORDS
INTERPRO
GOTERM_CC_DIRECT
GOTERM_BP_DIRECT
UP_KEYWORDS

Gap junction
IPR000990:Innexin
GO:0005921~gap junction
GO:0006811~ion transport
Cell junction

Annotation Cluster 60
Category
INTERPRO
GOTERM_MF_DIRECT
INTERPRO

Enrichment Score: 0.3254434616120558
Term
IPR016181:Acyl-CoA N-acyltransferase
GO:0008080~N-acetyltransferase activity
IPR000182:GNAT domain

Annotation Cluster 61
Category
COG_ONTOLOGY
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.31544957197360024
Term
Lipid metabolism / General function prediction only
IPR001206:Diacylglycerol kinase, catalytic domain
SM00046:DAGKc
IPR016064:ATP-NAD kinase-like domain

Annotation Cluster 62
Category
COG_ONTOLOGY
INTERPRO
INTERPRO

Enrichment Score: 0.31068677374002557
Term
Cell division and chromosome partitioning / Cytoskeleton
IPR000408:Regulator of chromosome condensation, RCC1

Annotation Cluster 63
Category
GOTERM_MF_DIRECT
INTERPRO
SMART

Enrichment Score: 0.30633464264590016
Term
GO:0005096~GTPase activator activity
IPR001164:Arf GTPase activating protein
SM00105:ArfGap

Annotation Cluster 64
Category
GOTERM_MF_DIRECT
INTERPRO
INTERPRO
UP_KEYWORDS
SMART
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 0.2891848053661439
Term
GO:0004930~G-protein coupled receptor activity
IPR017452:GPCR, rhodopsin-like, 7TM
IPR000276:G protein-coupled receptor, rhodopsin-like
G-protein coupled receptor
SM01381:SM01381
Receptor
Transducer

Annotation Cluster 65
Category
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.2786047593988897
Term
IPR013767:PAS fold
SM00091:PAS
IPR000014:PAS domain

Annotation Cluster 66
Category

Enrichment Score: 0.27226162298442
Term

155

INTERPRO
INTERPRO
INTERPRO

IPR001433:Oxidoreductase FAD/NAD(P)-binding
IPR017927:Ferredoxin reductase-type FAD-binding domain
IPR017938:Riboflavin synthase-like beta-barrel

Annotation Cluster 67
Category
UP_KEYWORDS
GOTERM_MF_DIRECT
INTERPRO

Enrichment Score: 0.27208541441613365
Term
Symport
GO:0005328~neurotransmitter:sodium symporter activity
IPR000175:Sodium:neurotransmitter symporter

Annotation Cluster 68
Category
GOTERM_BP_DIRECT
INTERPRO
INTERPRO

Enrichment Score: 0.2702830647661103
Term
GO:0008299~isoprenoid biosynthetic process
IPR000092:Polyprenyl synthetase
IPR008949:Terpenoid synthase

Annotation Cluster 69
Category
INTERPRO
INTERPRO
SMART
INTERPRO
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.2699629563427815
Term
IPR013098:Immunoglobulin I-set
IPR003598:Immunoglobulin subtype 2
SM00408:IGc2
IPR003599:Immunoglobulin subtype
IPR007110:Immunoglobulin-like domain
IPR013783:Immunoglobulin-like fold
SM00409:IG

Annotation Cluster 70
Category
INTERPRO
SMART
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.26550020870267177
Term
IPR020843:Polyketide synthase, enoylreductase
SM00829:SM00829
IPR002085:Alcohol dehydrogenase superfamily, zinc-type
IPR013154:Alcohol dehydrogenase GroES-like
IPR011032:GroES-like

Annotation Cluster 71
Category
INTERPRO
GOTERM_MF_DIRECT
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.2619741382632584
Term

Annotation Cluster 72
Category
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.251991370028117
Term

Annotation Cluster 73

Enrichment Score: 0.23337729281327935

GO:0004722~protein serine/threonine phosphatase activity
IPR015655:Protein phosphatase 2C
SM00332:PP2Cc
IPR001932:Protein phosphatase 2C (PP2C)-like

IPR002172:Low-density lipoprotein (LDL) receptor class A repeat
SM00192:LDLa

156

Category
INTERPRO
INTERPRO
INTERPRO

Term
IPR000873:AMP-dependent synthetase/ligase
IPR020845:AMP-binding, conserved site
IPR025110:Domain of unknown function DUF4009

Annotation Cluster 74
Category
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 0.23026181319429512
Term
Transit peptide
Mitochondrion
Mitochondrion inner membrane

Annotation Cluster 75
Category
GOTERM_MF_DIRECT
UP_KEYWORDS
INTERPRO
INTERPRO
INTERPRO
SMART
GOTERM_MF_DIRECT

Enrichment Score: 0.2271800792941584
Term

Annotation Cluster 76
Category
INTERPRO
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.22544319749190608
Term
IPR011765:Peptidase M16, N-terminal
IPR011237:Peptidase M16 domain
IPR007863:Peptidase M16, C-terminal domain
IPR011249:Metalloenzyme, LuxS/M16 peptidase-like

Annotation Cluster 77
Category
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.21006627649126675
Term
IPR008984:SMAD/FHA domain
IPR000253:Forkhead-associated (FHA) domain
SM00240:FHA

Annotation Cluster 78
Category
INTERPRO
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.20856794298872536
Term
IPR014710:RmlC-like jelly roll fold
IPR018490:Cyclic nucleotide-binding-like
IPR000595:Cyclic nucleotide-binding domain
SM00100:cNMP

Annotation Cluster 79
Category
INTERPRO
SMART
UP_KEYWORDS
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.19850784761592796
Term
IPR000033:LDLR class B repeat
SM00135:LY
EGF-like domain
IPR011042:Six-bladed beta-propeller, TolB-like
IPR000742:Epidermal growth factor-like domain
SM00181:EGF

Tyrosine-protein kinase
IPR020635:Tyrosine-protein kinase, catalytic domain
IPR008266:Tyrosine-protein kinase, active site
SM00219:TyrKc
GO:0004713~protein tyrosine kinase activity

157

INTERPRO
INTERPRO
INTERPRO
INTERPRO
SMART

IPR013032:EGF-like, conserved site
IPR000152:EGF-type aspartate/asparagine hydroxylation site
IPR018097:EGF-like calcium-binding, conserved site
IPR001881:EGF-like calcium-binding
SM00179:EGF_CA

Annotation Cluster 80
Category
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.1802901776532616
Term
IPR004087:K Homology domain
SM00322:KH
IPR004088:K Homology domain, type 1

Annotation Cluster 81
Category
INTERPRO
INTERPRO
GOTERM_MF_DIRECT
SMART

Enrichment Score: 0.1801460591858231
Term
IPR000504:RNA recognition motif domain
IPR012677:Nucleotide-binding, alpha-beta plait
GO:0000166~nucleotide binding
SM00360:RRM

Annotation Cluster 82
Category
INTERPRO
INTERPRO
INTERPRO
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.17775599841826317
Term
IPR020472:G-protein beta WD-40 repeat
IPR001680:WD40 repeat
IPR019775:WD40 repeat, conserved site
IPR015943:WD40/YVTN repeat-like-containing domain
SM00320:WD40
IPR017986:WD40-repeat-containing domain

Annotation Cluster 83
Category
INTERPRO
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.16148128454678556
Term
IPR020902:Actin/actin-like conserved site
IPR004001:Actin, conserved site
IPR004000:Actin-related protein
SM00268:ACTIN

Annotation Cluster 84
Category
INTERPRO
INTERPRO
INTERPRO
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.16081630065951844
Term
IPR015415:Vps4 oligomerisation, C-terminal
IPR003960:ATPase, AAA-type, conserved site
IPR003959:ATPase, AAA-type, core
IPR003439:ABC transporter-like
IPR003593:AAA+ ATPase domain
SM00382:AAA

Annotation Cluster 85
Category
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.15164473051291671
Term
IPR018980:FERM, C-terminal PH-like domain
SM01196:SM01196
IPR014352:FERM/acyl-CoA-binding protein, 3-helical bundle

158

INTERPRO
INTERPRO
INTERPRO
SMART
GOTERM_CC_DIRECT

IPR019748:FERM central domain
IPR000299:FERM domain
IPR019749:Band 4.1 domain
SM00295:B41
GO:0005856~cytoskeleton

Annotation Cluster 86
Category
INTERPRO
INTERPRO
INTERPRO
SMART
GOTERM_MF_DIRECT

Enrichment Score: 0.15163681821618055
Term
IPR013087:Zinc finger C2H2-type/integrase DNA-binding domain
IPR007087:Zinc finger, C2H2
IPR015880:Zinc finger, C2H2-like
SM00355:ZnF_C2H2
GO:0003676~nucleic acid binding

Annotation Cluster 87
Category
UP_KEYWORDS
GOTERM_MF_DIRECT
GOTERM_BP_DIRECT
GOTERM_BP_DIRECT
GOTERM_CC_DIRECT
GOTERM_CC_DIRECT
GOTERM_CC_DIRECT

Enrichment Score: 0.1447026494466192
Term
Initiation factor
GO:0003743~translation initiation factor activity
GO:0001731~formation of translation preinitiation complex
GO:0006446~regulation of translational initiation
GO:0033290~eukaryotic 48S preinitiation complex
GO:0005852~eukaryotic translation initiation factor 3 complex
GO:0016282~eukaryotic 43S preinitiation complex

Annotation Cluster 88
Category
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.13266945000734792
Term
IPR001005:SANT/Myb domain
SM00717:SANT
IPR017884:SANT domain

Annotation Cluster 89
Category
UP_KEYWORDS
UP_KEYWORDS
INTERPRO
GOTERM_MF_DIRECT
INTERPRO
INTERPRO
INTERPRO
INTERPRO
INTERPRO
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS
GOTERM_CC_DIRECT
UP_KEYWORDS
GOTERM_CC_DIRECT

Enrichment Score: 0.12602538486178053
Term
Cell junction
Ion transport
IPR018000:Neurotransmitter-gated ion-channel, conserved site
IPR002394:Nicotinic acetylcholine receptor
IPR006201:Neurotransmitter-gated ion-channel
IPR006202:Neurotransmitter-gated ion-channel ligand-binding
Ion channel
Cell membrane
Synapse
Ligand-gated ion channel
GO:0030054~cell junction
Postsynaptic cell membrane
GO:0045211~postsynaptic membrane

159

Annotation Cluster 90
Category
SMART
INTERPRO
INTERPRO

Enrichment Score: 0.11824316826224841
Term
SM00324:RhoGAP
IPR000198:Rho GTPase-activating protein domain
IPR008936:Rho GTPase activation protein

Annotation Cluster 91
Category
UP_KEYWORDS
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.1122557843932693
Term
ANK repeat
IPR002110:Ankyrin repeat
IPR020683:Ankyrin repeat-containing domain
SM00248:ANK

Annotation Cluster 92
Category
GOTERM_MF_DIRECT
INTERPRO
GOTERM_BP_DIRECT

Enrichment Score: 0.0989335301525666
Term
GO:0016831~carboxy-lyase activity
IPR002129:Pyridoxal phosphate-dependent decarboxylase
GO:0019752~carboxylic acid metabolic process

Annotation Cluster 93
Category
INTERPRO
INTERPRO
GOTERM_MF_DIRECT
INTERPRO

Enrichment Score: 0.09851299268707042
Term
IPR005829:Sugar transporter, conserved site
IPR005828:General substrate transporter
GO:0022891~substrate-specific transmembrane transporter activity
IPR003663:Sugar/inositol transporter

Annotation Cluster 94
Category
UP_KEYWORDS
UP_KEYWORDS
GOTERM_BP_DIRECT
UP_KEYWORDS

Enrichment Score: 0.09768413675411028
Term
Transcription regulation
Transcription
GO:0006351~transcription, DNA-templated
DNA-binding

Annotation Cluster 95
Category
UP_KEYWORDS
UP_KEYWORDS
INTERPRO
GOTERM_MF_DIRECT

Enrichment Score: 0.0837907843821477
Term
Flavoprotein
FAD

Annotation Cluster 96
Category
GOTERM_MF_DIRECT
INTERPRO
SMART
GOTERM_BP_DIRECT

Enrichment Score: 0.07081323934904664
Term
GO:0005089~Rho guanyl-nucleotide exchange factor activity
IPR000219:Dbl homology (DH) domain
SM00325:RhoGEF
GO:0035023~regulation of Rho protein signal transduction

Annotation Cluster 97
Category

Enrichment Score: 0.05047884032128006
Term

GO:0050660~flavin adenine dinucleotide binding

160

GOTERM_MF_DIRECT
INTERPRO
UP_KEYWORDS
INTERPRO
INTERPRO
INTERPRO
INTERPRO
SMART
SMART

GO:0004386~helicase activity

Annotation Cluster 98
Category
UP_KEYWORDS
UP_KEYWORDS
GOTERM_MF_DIRECT

Enrichment Score: 0.04918923169182432
Term
Isomerase
Rotamase
GO:0003755~peptidyl-prolyl cis-trans isomerase activity

Annotation Cluster 99
Category
COG_ONTOLOGY
INTERPRO
INTERPRO

Enrichment Score: 0.03762514320883115
Term
Lipid metabolism
IPR002018:Carboxylesterase, type B
IPR019819:Carboxylesterase type B, conserved site

Annotation Cluster 100
Category
INTERPRO
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.03326704008907252
Term
IPR019734:Tetratricopeptide repeat
IPR011990:Tetratricopeptide-like helical
SM00028:TPR
IPR013026:Tetratricopeptide repeat-containing domain

Annotation Cluster 101
Category
KEGG_PATHWAY
KEGG_PATHWAY
KEGG_PATHWAY
KEGG_PATHWAY
KEGG_PATHWAY

Enrichment Score: 0.01841829322131753
Term
ame03420:Nucleotide excision repair
ame03440:Homologous recombination
ame03410:Base excision repair
ame03030:DNA replication
ame03430:Mismatch repair

Annotation Cluster 102
Category
UP_KEYWORDS
GOTERM_MF_DIRECT
INTERPRO
INTERPRO
INTERPRO
UP_KEYWORDS
SMART

Enrichment Score: 0.003912689477190467
Term
Homeobox
GO:0043565~sequence-specific DNA binding
IPR001356:Homeodomain
IPR009057:Homeodomain-like
IPR017970:Homeobox, conserved site
DNA-binding
SM00389:HOX

Annotation Cluster 103
Category

Enrichment Score: 0.00310129416944087
Term

Helicase
IPR014014:RNA helicase, DEAD-box type, Q motif
IPR011545:DNA/RNA helicase, DEAD/DEAH box type, N-terminal
IPR001650:Helicase, C-terminal
IPR014001:Helicase, superfamily 1/2, ATP-binding domain
SM00490:HELICc
SM00487:DEXDc

161

UP_KEYWORDS
UP_KEYWORDS
GOTERM_CC_DIRECT

Cytoskeleton
Microtubule
GO:0005874~microtubule

Annotation Cluster 104
Category
GOTERM_BP_DIRECT
GOTERM_CC_DIRECT
GOTERM_MF_DIRECT

Enrichment Score: 1.6099659718154894E-9
Term
GO:0006412~translation
GO:0005840~ribosome
GO:0003735~structural constituent of ribosome

Up-Regulated
Annotation Cluster 1
Category
KEGG_PATHWAY
GOTERM_MF_DIRECT
GOTERM_BP_DIRECT
GOTERM_CC_DIRECT
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 22.67329895357984
Term
ame03010:Ribosome
GO:0003735~structural constituent of ribosome
GO:0006412~translation
GO:0005840~ribosome
Ribonucleoprotein
Ribosomal protein

Annotation Cluster 2
Category
SMART
INTERPRO
INTERPRO
GOTERM_CC_DIRECT
GOTERM_BP_DIRECT

Enrichment Score: 2.5402277919420313
Term
SM00651:Sm
IPR001163:Ribonucleoprotein LSM domain
IPR010920:Like-Sm (LSM) domain
GO:0005732~small nucleolar ribonucleoprotein complex
GO:0000398~mRNA splicing, via spliceosome

Annotation Cluster 3
Category
GOTERM_CC_DIRECT
INTERPRO
INTERPRO

Enrichment Score: 1.9332490052024476
Term
GO:0015934~large ribosomal subunit
IPR008991:Translation protein SH3-like domain
IPR014722:Ribosomal protein L2 domain 2

Annotation Cluster 4
Category
INTERPRO
INTERPRO
GOTERM_MF_DIRECT
SMART

Enrichment Score: 1.815880938858488
Term
IPR000504:RNA recognition motif domain
IPR012677:Nucleotide-binding, alpha-beta plait
GO:0000166~nucleotide binding
SM00360:RRM

Annotation Cluster 5
Category
GOTERM_CC_DIRECT
GOTERM_BP_DIRECT
INTERPRO
GOTERM_MF_DIRECT
SMART

Enrichment Score: 1.625730862268895
Term
GO:0005576~extracellular region
GO:0006030~chitin metabolic process
IPR002557:Chitin binding domain
GO:0008061~chitin binding
SM00494:ChtBD2

162

Annotation Cluster 6
Category
SMART
INTERPRO
UP_KEYWORDS
INTERPRO
SMART
INTERPRO
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 1.3835193027092916
Term
SM00409:IG
IPR003599:Immunoglobulin subtype
Immunoglobulin domain
IPR007110:Immunoglobulin-like domain
SM00408:IGc2
IPR003598:Immunoglobulin subtype 2
IPR013783:Immunoglobulin-like fold
IPR013106:Immunoglobulin V-set
IPR013098:Immunoglobulin I-set

Annotation Cluster 7
Category
KEGG_PATHWAY
GOTERM_MF_DIRECT
KEGG_PATHWAY
KEGG_PATHWAY

Enrichment Score: 1.265004431125057
Term
ame03020:RNA polymerase
GO:0003899~DNA-directed RNA polymerase activity
ame00230:Purine metabolism
ame00240:Pyrimidine metabolism

Annotation Cluster 8
Category
INTERPRO
INTERPRO
INTERPRO
INTERPRO
GOTERM_MF_DIRECT
INTERPRO
COG_ONTOLOGY

Enrichment Score: 1.131040798336076
Term
IPR024079:Metallopeptidase, catalytic domain
IPR018497:Peptidase M13, C-terminal domain
IPR008753:Peptidase M13, N-terminal domain
IPR000718:Peptidase M13
GO:0004222~metalloendopeptidase activity
IPR001590:Peptidase M12B, ADAM/reprolysin
Posttranslational modification, protein turnover, chaperones

Annotation Cluster 9
Category
SMART
INTERPRO
INTERPRO

Enrichment Score: 1.130658683742131
Term
SM00401:ZnF_GATA
IPR000679:Zinc finger, GATA-type
IPR013088:Zinc finger, NHR/GATA-type

Annotation Cluster 10
Category
UP_KEYWORDS
GOTERM_CC_DIRECT
UP_KEYWORDS

Enrichment Score: 1.0098475187450022
Term
Nucleus
GO:0005634~nucleus
DNA-binding

Annotation Cluster 11
Category
GOTERM_BP_DIRECT
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 0.9577376672440894
Term
GO:0006351~transcription, DNA-templated
Transcription
Transcription regulation

Annotation Cluster 12
Category

Enrichment Score: 0.9469867847459708
Term

163

INTERPRO
SMART
GOTERM_MF_DIRECT

IPR006170:Pheromone/general odorant binding protein
SM00708:PhBP
GO:0005549~odorant binding

Annotation Cluster 13
Category
INTERPRO
GOTERM_MF_DIRECT
KEGG_PATHWAY
UP_KEYWORDS
UP_KEYWORDS
INTERPRO

Enrichment Score: 0.9007408634727754
Term
IPR001353:Proteasome, subunit alpha/beta
GO:0004298~threonine-type endopeptidase activity
ame03050:Proteasome
Threonine protease
Proteasome
IPR016050:Proteasome, beta-type subunit, conserved site

INTERPRO
GOTERM_BP_DIRECT

IPR023333:Proteasome B-type subunit
GO:0051603~proteolysis involved in cellular protein catabolic process

GOTERM_CC_DIRECT

GO:0005839~proteasome core complex

Annotation Cluster 14
Category
INTERPRO
INTERPRO

Enrichment Score: 0.7858097993262126
Term
IPR013524:Runt domain
IPR000040:Acute myeloid leukemia 1 protein (AML1)/Runt

INTERPRO

IPR012346:p53/RUNT-type transcription factor, DNA-binding domain

INTERPRO

IPR008967:p53-like transcription factor, DNA-binding

Annotation Cluster 15
Category
INTERPRO

Enrichment Score: 0.6925895493040148
Term
IPR018525:Mini-chromosome maintenance, conserved site

SMART
INTERPRO

SM00350:MCM
IPR001208:Mini-chromosome maintenance, DNA-dependent ATPase

Annotation Cluster 16
Category
UP_KEYWORDS
GOTERM_CC_DIRECT
UP_KEYWORDS
SMART
INTERPRO
INTERPRO
GOTERM_BP_DIRECT
UP_KEYWORDS
GOTERM_BP_DIRECT
KEGG_PATHWAY

Enrichment Score: 0.6797839875585538
Term
Extracellular matrix
GO:0005578~proteinaceous extracellular matrix
Wnt signaling pathway
SM00097:WNT1
IPR018161:Wnt protein, conserved site
IPR005817:Wnt
GO:0016055~Wnt signaling pathway
Developmental protein
GO:0007275~multicellular organism development
ame04310:Wnt signaling pathway

Annotation Cluster 17
Category

Enrichment Score: 0.6410156227205825
Term

164

UP_KEYWORDS
GOTERM_CC_DIRECT
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS

Respiratory chain
GO:0070469~respiratory chain
Mitochondrion inner membrane
Mitochondrion
Electron transport

Annotation Cluster 18
Category
UP_KEYWORDS
INTERPRO

Enrichment Score: 0.6324483589697669
Term
Cell membrane
IPR006029:Neurotransmitter-gated ion-channel transmembrane domain

INTERPRO

IPR006202:Neurotransmitter-gated ion-channel ligand-binding

INTERPRO
GOTERM_MF_DIRECT

IPR006201:Neurotransmitter-gated ion-channel
GO:0005230~extracellular ligand-gated ion channel activity

INTERPRO

IPR018000:Neurotransmitter-gated ion-channel, conserved site

INTERPRO
UP_KEYWORDS
GOTERM_CC_DIRECT
GOTERM_CC_DIRECT
UP_KEYWORDS
UP_KEYWORDS
INTERPRO
GOTERM_MF_DIRECT

IPR006028:Gamma-aminobutyric acid A receptor
Synapse
GO:0045202~synapse
GO:0030054~cell junction
Cell junction
Postsynaptic cell membrane
IPR002394:Nicotinic acetylcholine receptor
GO:0004889~acetylcholine-activated cation-selective channel activity

GOTERM_CC_DIRECT
INTERPRO

GO:0045211~postsynaptic membrane
IPR027361:Nicotinic acetylcholine-gated receptor, transmembrane domain

UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS

Ion channel
Ligand-gated ion channel
Ion transport
Transport

Annotation Cluster 19
Category
UP_KEYWORDS
UP_KEYWORDS
UP_SEQ_FEATURE

Enrichment Score: 0.6132021357488939
Term
Cleavage on pair of basic residues
Amidation
signal peptide

Annotation Cluster 20
Category
INTERPRO
GOTERM_CC_DIRECT
GOTERM_BP_DIRECT

Enrichment Score: 0.608900338508913
Term
IPR009053:Prefoldin
GO:0016272~prefoldin complex
GO:0006457~protein folding

Annotation Cluster 21
Category
SMART

Enrichment Score: 0.5999526779710416
Term
SM01057:SM01057

165

INTERPRO
INTERPRO
KEGG_PATHWAY

IPR001148:Alpha carbonic anhydrase
IPR023561:Carbonic anhydrase, alpha-class
ame00910:Nitrogen metabolism

Annotation Cluster 22
Category
UP_KEYWORDS
UP_KEYWORDS
GOTERM_BP_DIRECT
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 0.591646731049165
Term
Antibiotic
Antimicrobial
GO:0042742~defense response to bacterium
Immunity
Innate immunity

Annotation Cluster 23
Category
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.5860972080757149
Term
IPR019956:Ubiquitin subgroup
IPR000626:Ubiquitin
SM00213:UBQ

Annotation Cluster 24
Category
INTERPRO
GOTERM_MF_DIRECT
SMART
INTERPRO
INTERPRO
INTERPRO
INTERPRO
GOTERM_BP_DIRECT
UP_KEYWORDS
SMART
INTERPRO
GOTERM_MF_DIRECT
GOTERM_CC_DIRECT
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 0.5579361534291644
Term
IPR017975:Tubulin, conserved site
GO:0005200~structural constituent of cytoskeleton
SM00864:SM00864
IPR023123:Tubulin, C-terminal
IPR000217:Tubulin
IPR008280:Tubulin/FtsZ, C-terminal
IPR003008:Tubulin/FtsZ, GTPase domain
GO:0007017~microtubule-based process
Cytoskeleton
SM00865:SM00865
IPR018316:Tubulin/FtsZ, 2-layer sandwich domain
GO:0003924~GTPase activity
GO:0005874~microtubule
Microtubule
GTP-binding

Annotation Cluster 25
Category
UP_KEYWORDS
INTERPRO
GOTERM_BP_DIRECT
UP_KEYWORDS

Enrichment Score: 0.5191879138867295
Term
Translocation
IPR004217:Tim10/DDP family zinc finger
GO:0015031~protein transport
Protein transport

Annotation Cluster 26
Category
UP_KEYWORDS
GOTERM_BP_DIRECT
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 0.5124139398620224
Term
DNA recombination
GO:0006310~DNA recombination
DNA repair
DNA damage

166

Annotation Cluster 27
Category
GOTERM_MF_DIRECT
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.4338065435913677
Term
GO:0005319~lipid transporter activity
IPR015816:Vitellinogen, beta-sheet N-terminal
IPR015819:Lipid transport protein, beta-sheet shell
IPR001747:Lipid transport protein, N-terminal

Annotation Cluster 28
Category
INTERPRO
INTERPRO
COG_ONTOLOGY
INTERPRO
GOTERM_MF_DIRECT

Enrichment Score: 0.42386467373610426
Term
IPR002018:Carboxylesterase, type B
IPR019826:Carboxylesterase type B, active site
Lipid metabolism
IPR019819:Carboxylesterase type B, conserved site
GO:0016787~hydrolase activity

Annotation Cluster 29
Category
INTERPRO
INTERPRO
SMART
INTERPRO
GOTERM_MF_DIRECT
SMART
SMART
INTERPRO

Enrichment Score: 0.400740212682717
Term
IPR002464:DNA/RNA helicase, ATP-dependent, DEAH-box type,
conserved site
IPR011709:Domain of unknown function DUF1605
SM00847:SM00847
IPR007502:Helicase-associated domain
GO:0008026~ATP-dependent helicase activity
SM00490:HELICc
SM00487:DEXDc
IPR014001:Helicase, superfamily 1/2, ATP-binding domain

INTERPRO
INTERPRO
INTERPRO

IPR001650:Helicase, C-terminal
IPR014014:RNA helicase, DEAD-box type, Q motif
IPR011545:DNA/RNA helicase, DEAD/DEAH box type, N-terminal

INTERPRO

IPR000629:RNA helicase, ATP-dependent, DEAD-box, conserved site

Annotation Cluster 30
Category
GOTERM_CC_DIRECT
GOTERM_BP_DIRECT
GOTERM_MF_DIRECT

Enrichment Score: 0.39178387696985084
Term
GO:0000276~mitochondrial proton-transporting ATP synthase complex,
coupling factor F(o)
GO:0015986~ATP synthesis coupled proton transport
GO:0015078~hydrogen ion transmembrane transporter activity

Annotation Cluster 31
Category
SMART
INTERPRO
UP_KEYWORDS
SMART
INTERPRO
INTERPRO

Enrichment Score: 0.38898139582192814
Term
SM00082:LRRCT
IPR000483:Cysteine-rich flanking region, C-terminal
Leucine-rich repeat
SM00369:LRR_TYP
IPR003591:Leucine-rich repeat, typical subtype
IPR001611:Leucine-rich repeat

167

Annotation Cluster 32
Category
INTERPRO

Enrichment Score: 0.3648985169547037
Term
IPR009000:Translation elongation/initiation factor/Ribosomal, beta-barrel

INTERPRO

IPR004161:Translation elongation factor EFTu/EF1A, domain 2

GOTERM_MF_DIRECT
INTERPRO

GO:0003924~GTPase activity
IPR000795:Elongation factor, GTP-binding domain

Annotation Cluster 33
Category
INTERPRO
INTERPRO
GOTERM_MF_DIRECT
INTERPRO

Enrichment Score: 0.34621265468482626
Term
IPR011645:Haem NO binding associated
IPR001054:Adenylyl cyclase class-3/4/guanylyl cyclase
GO:0004383~guanylate cyclase activity
IPR018297:Adenylyl cyclase class-3/4/guanylyl cyclase, conserved site

SMART
UP_KEYWORDS
GOTERM_BP_DIRECT

SM00044:CYCc
Lyase
GO:0035556~intracellular signal transduction

Annotation Cluster 34
Category
SMART
INTERPRO
UP_KEYWORDS
SMART
INTERPRO

Enrichment Score: 0.3440696455289717
Term
SM00072:GuKc
IPR008145:Guanylate kinase/L-type calcium channel
SH3 domain
SM00326:SH3
IPR001452:Src homology-3 domain

Annotation Cluster 35
Category
INTERPRO

Enrichment Score: 0.3400209199980145
Term
IPR002172:Low-density lipoprotein (LDL) receptor class A repeat

SMART
INTERPRO

SM00192:LDLa
IPR023415:Low-density lipoprotein (LDL) receptor class A, conserved site

Annotation Cluster 36
Category
UP_KEYWORDS
GOTERM_MF_DIRECT
UP_KEYWORDS

Enrichment Score: 0.329492518239019
Term
Metalloprotease
GO:0008237~metallopeptidase activity
Glycoprotein

Annotation Cluster 37
Category
UP_KEYWORDS
INTERPRO
INTERPRO
SMART
UP_KEYWORDS

Enrichment Score: 0.32693725189266776
Term
G-protein coupled receptor
IPR017452:GPCR, rhodopsin-like, 7TM
IPR000276:G protein-coupled receptor, rhodopsin-like
SM01381:SM01381
Receptor

168

UP_KEYWORDS
GOTERM_MF_DIRECT

Transducer
GO:0004930~G-protein coupled receptor activity

Annotation Cluster 38
Category
INTERPRO
INTERPRO
UP_KEYWORDS
SMART
INTERPRO
GOTERM_MF_DIRECT
UP_KEYWORDS
INTERPRO
GOTERM_BP_DIRECT

Enrichment Score: 0.31768016228988083
Term
IPR020479:Homeodomain, metazoa
IPR009057:Homeodomain-like
DNA-binding
SM00389:HOX
IPR001356:Homeodomain
GO:0043565~sequence-specific DNA binding
Homeobox
IPR017970:Homeobox, conserved site
GO:0006355~regulation of transcription, DNA-templated

Annotation Cluster 39
Category
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.3172230438459506
Term
IPR001304:C-type lectin
IPR016186:C-type lectin-like
IPR016187:C-type lectin fold

Annotation Cluster 40
Category
SMART
INTERPRO
INTERPRO

Enrichment Score: 0.30358499885476903
Term
SM00339:FH
IPR001766:Transcription factor, fork head
IPR011991:Winged helix-turn-helix DNA-binding domain

Annotation Cluster 41
Category
SMART
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.30133995954774584
Term
SM00449:SPRY
IPR003877:SPla/RYanodine receptor SPRY
IPR001870:B30.2/SPRY domain
IPR013320:Concanavalin A-like lectin/glucanase, subgroup

Annotation Cluster 42
Category
INTERPRO
SMART
INTERPRO
INTERPRO

Enrichment Score: 0.28921563699416075
Term
IPR001314:Peptidase S1A, chymotrypsin-type
SM00020:Tryp_SPc
IPR001254:Peptidase S1
IPR009003:Trypsin-like cysteine/serine peptidase domain

GOTERM_MF_DIRECT
INTERPRO
UP_KEYWORDS

GO:0004252~serine-type endopeptidase activity
IPR018114:Peptidase S1, trypsin family, active site
Serine protease

Annotation Cluster 43
Category

Enrichment Score: 0.2527851472350029
Term

169

SMART
INTERPRO
INTERPRO

SM00298:CHROMO
IPR016197:Chromo domain-like
IPR000953:Chromo domain/shadow

Annotation Cluster 44
Category
INTERPRO
GOTERM_CC_DIRECT
UP_KEYWORDS
UP_KEYWORDS
INTERPRO

Enrichment Score: 0.24026228886816542
Term
IPR009072:Histone-fold
GO:0000786~nucleosome
Chromosome
Nucleosome core
IPR007125:Histone core

Annotation Cluster 45
Category
INTERPRO
INTERPRO
INTERPRO

KEGG_PATHWAY
GOTERM_BP_DIRECT
KEGG_PATHWAY
GOTERM_BP_DIRECT

Enrichment Score: 0.23476796871215266
Term
IPR001557:L-lactate/malate dehydrogenase
IPR022383:Lactate/malate dehydrogenase, C-terminal
IPR015955:Lactate dehydrogenase/glycoside hydrolase, family 4, Cterminal
GO:0016616~oxidoreductase activity, acting on the CH-OH group of
donors, NAD or NADP as acceptor
ame00270:Cysteine and methionine metabolism
GO:0019752~carboxylic acid metabolic process
ame00620:Pyruvate metabolism
GO:0005975~carbohydrate metabolic process

Annotation Cluster 46
Category
GOTERM_MF_DIRECT
INTERPRO
INTERPRO

Enrichment Score: 0.22385608725379771
Term
GO:0046872~metal ion binding
IPR007087:Zinc finger, C2H2
IPR013087:Zinc finger C2H2-type/integrase DNA-binding domain

SMART
INTERPRO

SM00355:ZnF_C2H2
IPR015880:Zinc finger, C2H2-like

Annotation Cluster 47
Category
INTERPRO
INTERPRO
GOTERM_BP_DIRECT

Enrichment Score: 0.21992443018159333
Term
IPR000488:Death domain
IPR011029:Death-like domain
GO:0007165~signal transduction

Annotation Cluster 48
Category
SMART
INTERPRO
INTERPRO
INTERPRO
INTERPRO

Enrichment Score: 0.21316119066385328
Term
SM00320:WD40
IPR001680:WD40 repeat
IPR019775:WD40 repeat, conserved site
IPR017986:WD40-repeat-containing domain
IPR015943:WD40/YVTN repeat-like-containing domain

Annotation Cluster 49
Category

Enrichment Score: 0.20816885012417152
Term

GOTERM_MF_DIRECT

170

INTERPRO
PIR_SUPERFAMILY
INTERPRO

IPR012132:Glucose-methanol-choline oxidoreductase
PIRSF000137:glucose-methanol-choline oxidoreductase
IPR007867:Glucose-methanol-choline oxidoreductase, C-terminal

INTERPRO

IPR000172:Glucose-methanol-choline oxidoreductase, N-terminal

GOTERM_MF_DIRECT

GO:0016614~oxidoreductase activity, acting on CH-OH group of donors

GOTERM_MF_DIRECT
INTERPRO

GO:0050660~flavin adenine dinucleotide binding
IPR023753:Pyridine nucleotide-disulphide oxidoreductase, FAD/NAD(P)binding domain

Annotation Cluster 50
Category
INTERPRO
GOTERM_MF_DIRECT
UP_KEYWORDS
INTERPRO

Enrichment Score: 0.1791919509442036
Term
IPR005814:Aminotransferase class-III
GO:0008483~transaminase activity
Pyridoxal phosphate
IPR015421:Pyridoxal phosphate-dependent transferase, major region,
subdomain 1
GO:0030170~pyridoxal phosphate binding
IPR015424:Pyridoxal phosphate-dependent transferase
IPR015422:Pyridoxal phosphate-dependent transferase, major region,
subdomain 2

GOTERM_MF_DIRECT
INTERPRO
INTERPRO

Annotation Cluster 51
Category
SMART
INTERPRO

Enrichment Score: 0.17779540327888022
Term
SM01100:SM01100
IPR001071:Cellular retinaldehyde binding/alpha-tocopherol transport

INTERPRO
SMART
INTERPRO
GOTERM_MF_DIRECT

IPR011074:CRAL/TRIO, N-terminal domain
SM00516:SEC14
IPR001251:CRAL-TRIO domain
GO:0005215~transporter activity

Annotation Cluster 52
Category
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS

Enrichment Score: 0.16468810755496266
Term
Zinc
Metal-binding
Zinc-finger

Annotation Cluster 53
Category
GOTERM_MF_DIRECT
UP_KEYWORDS
SMART
INTERPRO
INTERPRO
INTERPRO
UP_KEYWORDS
GOTERM_MF_DIRECT
INTERPRO

Enrichment Score: 0.15292023006330072
Term
GO:0004672~protein kinase activity
Serine/threonine-protein kinase
SM00220:S_TKc
IPR000719:Protein kinase, catalytic domain
IPR008271:Serine/threonine-protein kinase, active site
IPR017441:Protein kinase, ATP binding site
Kinase
GO:0004674~protein serine/threonine kinase activity
IPR011009:Protein kinase-like domain

171

UP_KEYWORDS

Transferase

Annotation Cluster 54
Category
INTERPRO

Enrichment Score: 0.14849403552329868
Term
IPR019594:Glutamate receptor, L-glutamate/glycine-binding

GOTERM_MF_DIRECT

GO:0005234~extracellular-glutamate-gated ion channel activity

INTERPRO
GOTERM_MF_DIRECT

IPR001320:Ionotropic glutamate receptor
GO:0004970~ionotropic glutamate receptor activity

Annotation Cluster 55
Category
GOTERM_MF_DIRECT
INTERPRO
GOTERM_BP_DIRECT
GOTERM_CC_DIRECT
INTERPRO

Enrichment Score: 0.1465758655346196
Term
GO:0015035~protein disulfide oxidoreductase activity
IPR012336:Thioredoxin-like fold
GO:0045454~cell redox homeostasis
GO:0005623~cell
IPR013766:Thioredoxin domain

Annotation Cluster 56
Category
INTERPRO
SMART
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.14457633756131663
Term
IPR000742:Epidermal growth factor-like domain
SM00179:EGF_CA
IPR001881:EGF-like calcium-binding
SM00181:EGF
IPR000152:EGF-type aspartate/asparagine hydroxylation site

INTERPRO
INTERPRO
UP_KEYWORDS
INTERPRO

IPR018097:EGF-like calcium-binding, conserved site
IPR013032:EGF-like, conserved site
EGF-like domain
IPR009030:Insulin-like growth factor binding protein, N-terminal

GOTERM_MF_DIRECT

GO:0005509~calcium ion binding

Annotation Cluster 57
Category
GOTERM_MF_DIRECT
SMART
INTERPRO
INTERPRO
UP_KEYWORDS

Enrichment Score: 0.1331274369109913
Term
GO:0005216~ion channel activity
SM00248:ANK
IPR002110:Ankyrin repeat
IPR020683:Ankyrin repeat-containing domain
ANK repeat

Annotation Cluster 58
Category
UP_KEYWORDS
INTERPRO

Enrichment Score: 0.13063392323374434
Term
Isomerase
IPR002130:Cyclophilin-like peptidyl-prolyl cis-trans isomerase domain

INTERPRO

IPR024936:Cyclophilin-type peptidyl-prolyl cis-trans isomerase

UP_KEYWORDS

Rotamase

172

GOTERM_BP_DIRECT
GOTERM_MF_DIRECT

GO:0006457~protein folding
GO:0003755~peptidyl-prolyl cis-trans isomerase activity

Annotation Cluster 59
Category
GOTERM_MF_DIRECT
INTERPRO
INTERPRO

Enrichment Score: 0.12260609909013652
Term
GO:0008080~N-acetyltransferase activity
IPR000182:GNAT domain
IPR016181:Acyl-CoA N-acyltransferase

Annotation Cluster 60
Category
SMART
INTERPRO
INTERPRO

Enrichment Score: 0.11097685475483653
Term
SM00322:KH
IPR004087:K Homology domain
IPR004088:K Homology domain, type 1

Annotation Cluster 61
Category
INTERPRO
SMART
UP_KEYWORDS

Enrichment Score: 0.09367797840987173
Term
IPR001781:Zinc finger, LIM-type
SM00132:LIM
LIM domain

Annotation Cluster 62
Category
UP_KEYWORDS
GOTERM_MF_DIRECT
COG_ONTOLOGY

Enrichment Score: 0.08458418914700383
Term
Monooxygenase
GO:0020037~heme binding
Secondary metabolites biosynthesis, transport, and catabolism

INTERPRO
UP_KEYWORDS
INTERPRO
GOTERM_MF_DIRECT

IPR002401:Cytochrome P450, E-class, group I
Heme
IPR001128:Cytochrome P450
GO:0016705~oxidoreductase activity, acting on paired donors, with
incorporation or reduction of molecular oxygen

INTERPRO
GOTERM_MF_DIRECT
UP_KEYWORDS
GOTERM_MF_DIRECT
UP_KEYWORDS

IPR017972:Cytochrome P450, conserved site
GO:0004497~monooxygenase activity
Iron
GO:0005506~iron ion binding
Oxidoreductase

Annotation Cluster 63
Category
SMART
INTERPRO
INTERPRO

Enrichment Score: 0.08261573613344135
Term
SM00454:SAM
IPR001660:Sterile alpha motif domain
IPR013761:Sterile alpha motif/pointed domain

Annotation Cluster 64
Category
INTERPRO
INTERPRO

Enrichment Score: 0.08249855974494377
Term
IPR005828:General substrate transporter
IPR003663:Sugar/inositol transporter

173

GOTERM_MF_DIRECT

GO:0022891~substrate-specific transmembrane transporter activity

INTERPRO

IPR005829:Sugar transporter, conserved site

Annotation Cluster 65
Category
UP_KEYWORDS
INTERPRO
GOTERM_CC_DIRECT
GOTERM_MF_DIRECT
INTERPRO
GOTERM_BP_DIRECT
INTERPRO

Enrichment Score: 0.05823738247560854
Term
GTP-binding
IPR020849:Small GTPase superfamily, Ras type
GO:0016020~membrane
GO:0005525~GTP binding
IPR005225:Small GTP-binding protein domain
GO:0007264~small GTPase mediated signal transduction
IPR001806:Small GTPase superfamily

Annotation Cluster 66
Category
INTERPRO
INTERPRO
GOTERM_BP_DIRECT

Enrichment Score: 0.04909501010032386
Term
IPR011701:Major facilitator superfamily
IPR020846:Major facilitator superfamily domain
GO:0055085~transmembrane transport

Annotation Cluster 67
Category
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.03187953124700801
Term
IPR000210:BTB/POZ-like
IPR011333:BTB/POZ fold
SM00225:BTB

Annotation Cluster 68
Category
INTERPRO
INTERPRO
INTERPRO
SMART

Enrichment Score: 0.026813697036988098
Term
IPR003960:ATPase, AAA-type, conserved site
IPR003959:ATPase, AAA-type, core
IPR003593:AAA+ ATPase domain
SM00382:AAA

Annotation Cluster 69
Category
UP_KEYWORDS
UP_KEYWORDS
GOTERM_MF_DIRECT

Enrichment Score: 0.01851786526800038
Term
Nucleotide-binding
ATP-binding
GO:0005524~ATP binding

Annotation Cluster 70
Category
GOTERM_CC_DIRECT
GOTERM_MF_DIRECT
UP_KEYWORDS
UP_KEYWORDS
INTERPRO
GOTERM_MF_DIRECT

Enrichment Score: 0.016374083153132487
Term
GO:0005886~plasma membrane
GO:0005549~odorant binding
Olfaction
Sensory transduction
IPR004117:Olfactory receptor, Drosophila
GO:0004984~olfactory receptor activity

Annotation Cluster 71
Category

Enrichment Score: 0.009703276678371543
Term

174

INTERPRO
INTERPRO
GOTERM_MF_DIRECT
INTERPRO

IPR018247:EF-Hand 1, calcium-binding site
IPR011992:EF-hand-like domain
GO:0005509~calcium ion binding
IPR002048:EF-hand domain

Annotation Cluster 72
Category
INTERPRO
INTERPRO
SMART
INTERPRO

Enrichment Score: 0.004441392633435781
Term
IPR013026:Tetratricopeptide repeat-containing domain
IPR011990:Tetratricopeptide-like helical
SM00028:TPR
IPR019734:Tetratricopeptide repeat

Annotation Cluster 73
Category
UP_KEYWORDS
UP_KEYWORDS
UP_KEYWORDS
GOTERM_CC_DIRECT

Enrichment Score: 8.697889380509524E-5
Term
Membrane
Transmembrane helix
Transmembrane
GO:0016021~integral component of membrane

175

Appendix S: Enriched Genes at Control-4h and Treatment-0h. DAVID was used for
enrichment analysis. Genes that had an FDR<0.05 were considered significantly enriched. Downregulated genes had an FDR <0.09. Included in this table for comparison.
Annotation
Number

Function of
Cluster

Enrichment
Score

1

Transmembrane

0.7

1

Transmembrane

4.59

Upregulated

√

176

DownRegulated

Unique Genes in Cluster

√

uncharacterized protein
PF11_0207-like
(LOC100578519)
uncharacterized LOC725682
sex peptide receptor-like
(LOC724225)
dual oxidase (LOC551970)
beta-1,4-glucuronyltransferase 1like(LOC727358)
CAAX prenyl protease 1
homolog(LOC551466)
cadherin-23(LOC410368)
dopamine receptor, D1(Dop1)
facilitated trehalose transporter
Tret1-like(LOC100577385)
glutamate receptor 3like(LOC100577496)
GPI mannosyltransferase
4(LOC551642)
innexin inx7(LOC552285)
leukocyte tyrosine kinase
receptor(LOC408718)
nanchung(LOC552792)
neprilysin-2(LOC724616)
neuropeptide Y receptorlike(LOC411614)
odorant receptor 27(Or27)
odorant receptor 30(Or30)
organic cation transporter
protein-like(LOC102653922)
probable phospholipidtransporting ATPase
VD(LOC409192)
proton-coupled amino acid
transporter 4(LOC410741)
sodium/bile acid cotransporter 7like(LOC724999)
sushi, von Willebrand factor type
A, EGF and pentraxin domaincontaining
protein 1like(LOC100578434)
transmembrane protein
114(LOC725318)
uncharacterized
LOC100576461(LOC100576461)
uncharacterized
LOC100576683(LOC100576683)

2

Immunoglobulin

1.15

√

177

uncharacterized
LOC551778(LOC551778)
uncharacterized
LOC725724(LOC725724)
uncharacterized
LOC727444(LOC727444)
zinc transporter 1(LOC727479)
uncharacterized LOC100576471
uncharacterized LOC725354
cell adhesion molecule 2-like
(LOC409000)

178