Pesticide and pathogen exposure causes idiosyncratic gene expression responses across four diverse North American bumble bee species

Bumble bee populations of certain species have been declining precipitously in North America over several decades. Hypotheses for declines include exposure to the fungal pathogen Nosema bombi and neonicotinoid pesticides. Importantly, populations of some bumble bee species remain stable despite their presumed exposure to these same stressors. Here we hypothesize that declining and stable species, respectively, exhibit distinct responses to N. bombi and neonicotinoid pesticides, detectable as differential gene expression profiles. To test this, we exposed developing larvae from colonies of Bombus occidentalis (declining) and B. impatiens (stable) to N. bombi and to the neonicotinoid imidacloprid, plus a combination of both. RNA-seq analysis revealed almost no overlap between these species in gene expression responses to the individual stressors. Combined-stressor effects were more similar, but nevertheless differed significantly in differentially expressed genes and in gene coexpression network modules. To test whether the molecular responses correlated within declining and stable species, we carried out quantitative PCR on twenty target genes and included a second pair of species, B. terricola (declining) and B. griseocollis (stable). Results showed no correlation with decline or stable status, with each of the four species exhibiting species-specific responses. Overall, these results highlight that generalizing causes of decline across different species may be misleading, as diverse species respond in a species-specific manner to particular environmental stressors.

Methods
Experimental design
Colonies from four North American bumble bee species were set in the Illinois State University facilities: Bombus impatiens (n = 7), Bombus occidentalis (n = 3), Bombus griseocollis (n = 4) and Bombus terricola (n = 2). Four microcolonies consisting of five adult workers and a larval brood clump were set from each parental colony. Each set of microcolonies from a source colony was exposed to one of four possible treatments:

unexposed to imidacloprid and N. bombi (control)
exposed to N. bombi, unexposed to imidacloprid
exposed to imidacloprid, unexposed to N. bombi

exposed to both imidacloprid and N. bombi

After 24 h of provisioning the microcolonies with untreated sugar water and pollen ad libitum, both were replaced with either control or imidacloprid (7 ppb) treated provisions for 72 h. After the first exposure event, provisions were again switched to untreated, and larvae were individually fed a 2 μL inoculum of either control or Nosema bombi (20,000 spores). After 24 h of the second exposure event, four larvae from each microcolony was flash frozen in liquid nitrogen and stored at -80 ºC.

RNA extraction
Total RNA was extracted from three larvae (leaving the fourth as backup) from each set of microcolonies with the E.Z.N.A. Total RNA Kit I (Omega Bio-tek), followoning manufacturers indications and adding a DNAse I treatment. Total RNA integrity was checked by running 1 μL in a 1% agarose gel, and quantity was assessed by analyzing 1 μL with a Qubit fluorometer.

RNA sequencing and bioinformatic analysis
For B. impatiens and B. occidentalis samples, RNA from the three larvae per microcolony was pooled in a single tube to obtain 1 μg of total RNA. Library contruction involved poly-A tail selection, followed by sequencing by a Illumina NovaSeq6000 sequencer. Single end raw reads are available in the SRA repository (NCBI) under the BioProject PRJNA1071869.

Raw reads were trimmed and cleaned from adapter sequences with Trimmomatic v0.38. FastQC was used to check read quality before and after the trimming. Trimmed reads were mapped to the B. impatiens genome v2.2 (GenBank accession number GCF_000188095.3) with STAR v2.7. Read counts were summarized per gene feature with htseq-count. Differential gene expression analysis was performed with DESeq2. Co-expression networks were built with WGCNA. Gene ontology term enrichment analysis was performed with topGO. The RNA-seq analysis pipeline and instructions on how to run it are described in the README file.

Quantitative PCR analysis
Primers for 20 differentially expressed genes were designed with NCBI's PrimerBLAST, using the available bumble bee genomes as templates, to test in all four bumble bee species. From 1 µg of RNA per sample, cDNA was synthesized with the iScript cDNA kit (BioRad). Quantitative qPCR was performed with the Luna qPCR reaction kit (NEBL) using 25 ng of cDNA per reaction as template. Reactions were run in a StepOne Thermocycler (Thermo Fisher Scientific), with the following conditions: hot start at 95 °C for 1 minute, then 40 cycles at 95 °C for 15 seconds and 60 °C for 30 seconds each, and finally a melt curve ranging from 60 °C to 95 °C with reads occurring every 0.3 °C. Stability of housekeeping genes was analyzed with GeNorm v3. Relative quantities were calculated following the Pfaffl method. Details on running the qPCR analysis are described in the README file.