# RNA-sequence analysis of gene expression from honeybees (Apis mellifera) infected with Nosema ceranae

Bouabid Badaoui, Andre Fougeroux, Fabien Petit, Anna Anselmo, Chiara Gorni, Marco Cucurachi, Antonella Cersini, Anna Granato, Giusy Cardeti, Giovanni Formato, et al.

► To cite this version:

Bouabid Badaoui, Andre Fougeroux, Fabien Petit, Anna Anselmo, Chiara Gorni, et al. RNA-sequence analysis of gene expression from honeybees (Apis mellifera) infected with Nosema ceranae. PLoS ONE, 2017, 12 (3), 10.1371/journal.pone.0173438. [DOI: 10.1371/journal.pone.0173438](https://doi.org/10.1371/journal.pone.0173438)

HAL Id: hal-01604257

https://hal.science/hal-01604257

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# RNA-sequence analysis of gene expression from honeybees (*Apis mellifera*) infected with *Nosema ceranae*

Bouabid Badaoui¹ na, André Fougeroux², Fabien Petit³, Anna Anselmo¹, Chiara Gorni¹, Marco Cucurachi¹, Antonella Cersini⁴, Anna Granato⁵, Giusy Cardeti⁴, Giovanni Formato⁴, Franco Mutinelli⁵, Elisabetta Giuffra¹ nb, John L. Williams⁶, Sara Botti¹*

1 Parco Tecnologico Padano - CERSA, Integrative Biology Group, Lodi, Italy, 2 Syngenta, Guyancourt, France, 3 RNAGRO, TOURS, France, 4 Istituto Zooprofilattico Sperimentale del Lazio e della Toscana, Roma, Italy, 5 Istituto Zooprofilattico Sperimentale delle Venezie, Legnaro, Padua, Italy, 6 Davies Research Centre, University of Adelaide, Roseworthy, South Australia, Australia

na Current address: Mohamed V University in Rabat, Science Faculty, Rabat, Morocco

nb Current address: INRA - laboratory of Animal Genetics and Integrative Biology Team Genetics, Domaine de Vilvert, Jouy-en-Josas, France

* sara.botti@ptp.it

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Badaoui B, Fougeroux A, Petit F, Anselmo A, Gorni C, Cucurachi M, et al. (2017) RNA-sequence analysis of gene expression from honeybees (*Apis mellifera*) infected with *Nosema ceranae*. PLoS ONE 12(3): e0173438. https://doi.org/10.1371/journal.pone.0173438

**Editor:** Olav Rueppell, University of North Carolina at Greensboro, UNITED STATES

**Received:** September 2, 2015

**Accepted:** February 22, 2017

**Published:** March 28, 2017

**Copyright:** © 2017 Badaoui et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** The data may be found at the SRA database at the following URL: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA378655/ and accession number PRJNA378655.

**Funding:** Funded by Eureka-Eurostars Project Application E!5069 NewBone. Andrea Fougeroux was employed by Syngenta during the course of the study. Syngenta provided support in the form of salary for author AF, but did not have any additional role in the study design, data collection

## Abstract

Honeybees (*Apis mellifera*) are constantly subjected to many biotic stressors including parasites. This study examined honeybees infected with *Nosema ceranae* (*N. ceranae*). *N. ceranae* infection increases the bees energy requirements and may contribute to their decreased survival. RNA-seq was used to investigate gene expression at days 5, 10 and 15 Post Infection (P.I) with *N. ceranae*. The expression levels of genes, isoforms, alternative transcription start sites (TSS) and differential promoter usage revealed a complex pattern of transcriptional and post-transcriptional gene regulation suggesting that bees use a range of tactics to cope with the stress of *N. ceranae* infection. *N. ceranae* infection may cause reduced immune function in the bees by: (i) disturbing the host amino acids metabolism (ii) down-regulating expression of antimicrobial peptides (iii) down-regulation of cuticle coatings and (iv) down-regulation of odorant binding proteins.

## Introduction

Honeybees (*Apis mellifera*), are critical for agricultural ecosystems, but are exposed to many biotic stressors including bacteria, viruses and fungi. *Nosema ceranae* (*N. ceranae*) is a microsporidian parasite which is one of the most common parasites of the honeybee. It infects the midgut [1] and has many adverse effects which reduce the health of individual honeybees [2–4] and of the colony [5–10]. *N. ceranae* infection is associated with several fitness related problems including deterioration of the gut epithelial cells [11,12], immune suppression [13] and energetic stress [14,15]. Honeybees have developed many defense mechanisms to deal with parasite infection which include an antimicrobial secretion layer on the gut, and cellular and
and analysis, decision to publish, or preparation of
the manuscript. The specific role of this author is
articulated in the 'author contributions' section.

**Competing interests:** Andrea Fougeroux was employed by Syngenta during the course of the study. Fabien Petit acted as a contractor to Syngenta during the course of the study. There are no patents, products in development or marketed products to declare. This does not alter our adherence to all the PLOS ONE policies on sharing data and materials.

humoral immune mechanisms [16]. The honeybee humoral response has been shown to involve four antimicrobial peptides, abaecin, apidaecin, defensin and hymenoptaecin [17–20].

A draft *A. mellifera* genome sequence was published in 2004 [21] and an improved version of the genome sequence was published in 2014 [22]. This sequence is valuable for the interpretation of RNA-seq data to study gene expression. RNA-seq may be used [23] to characterize the whole transcriptome including transcription start sites (TSS), splicing variants and differential promoter usage [24].

The molecular response of the honey bee to infection by *N. ceranae* has been investigated using RNA-seq, but the focus to date has been on cells of the midgut [25]. Digital gene expression has also been used to study the neurogenomic response to *N. ceranae* in the bee brain [26]. The honey-bee nervous system, trachea and intestine interact to maintain the physiological homeostasis following damage resulting from infection; e.g. intestinal stem cells replace damaged intestinal epithelial cells [27]. Interestingly, the message for intestinal regeneration come not only from the damaged epithelial cells, but also from adjacent (trachea, muscles) and distant tissues, and the brain (reviewed in [28]).

In this study, RNA-seq was used to investigate gene expression, splice variants, TSSs and differential promoter usage in honeybees infected with *N. ceranae* at days 5, 10 and 15 P.I. The new version of the bee genome (Amel_4.5_scaffolds.fa) was used to annotate the genes. The study addressed the global transcriptome response of the bee to *N. ceranae* infection, not the response of specific tissues, and identified many transcripts and pathways involved in the interaction between the honeybee and the parasite.

**Material and methods**

*N. ceranae* spore preparation

*N. ceranae* spores were obtained according to Aufauvre et al. 2011 [25]. Briefly, forager bees were sampled at the hive entrance. The abdomens of 30 bees were triturated adding 30 ml distilled water. The homogenate was filtered and the resulting suspension was centrifuged for 6 minutes at 800g. The supernatant was removed and the spores resuspended in 3 ml distilled water. The spore concentration was determined by counting using a haemocytometer chamber. A median of 8.6 X 10⁵ spores/ml was obtained.

Inoculation of bees with *N. ceranae* and RNA extraction from honeybees

Frames of *Apis mellifera* capped brood containing "about-to-emerge" pupae were cut out from healthy colonies located in France. The hives were verified as *N. ceranae* free. Subsequently, the frame of pupae from each colony was divided into 2 equal parts, introduced into two plastic boxes and kept in an incubator at ~30°C with high humidity. The experiment was performed in the apiary of Luz Saint Sauveur located at 00° 00' 50.4"O and 42° 51' 57.6" N.

After about 2 days, the emerging bees were counted, and 30 bees were retained within each box. The two boxes were transferred to the inoculation chamber (25°C) with a 12 hour dark/light (12/12 d/l) cycle. The bees of the "inoculation box" were fed for two days with 2 ml of 66% sucrose solution (w/v) containing 86,000 *N. ceranae* spores per bee (total of 2.6X10⁵ spores per "inoculation box"). The bees of the "control box" were fed with 2 ml of a solution of 66% sucrose without *N. ceranae* spores for two days. Throughout the experiment, sucrose solution was supplemented only when the Pasteur pipettes in the boxes were completely empty. On days 5, 10 and 15, five bees per box were removed and ground in liquid nitrogen, transferred in micro-tube containing Trizol and stored at -80°C. The number of live/dead bees in both boxes was counted daily and the dead bees were carefully removed. Total RNA was extracted from three bees from each box and at each time point (5, 10 and 15 days post
infection, P.I.) using TRIzol (Invitrogen Life Technologies, Milan-Italy) and RNeasy columns (Qiagen). RNA quality was assessed by microcapillary electrophoresis on an Agilent 2001 Bioanalyzer (Agilent Technologies) with RNA 6000 Nanochips. RNA was quantified by spectrophotometry (ND-1000; NanoDrop Technologies).

**Control for N. ceranae infection.** Presence of *N. ceranae* RNA in each sample was assessed by PCR [2] or qPCR [4,29]. Only samples from the "inoculation box" that were positive for *N. ceranae* RNA were used to create RNA sequencing libraries. To verify that the inoculation was successful, the 15 remaining bees after 14 days (P.I) were examined for *N. ceranae* spores by microscopy and the number of spores present recorded.

## Library preparation for RNA-seq

Eighteen libraries corresponding to 9 control and 9 infected samples at 5, 10 and 15 days P.I. (3 samples each time-point) were sequenced. TruSeq Sample Prep Kits (Illumina Inc., San Diego, CA) were used to create libraries for 1×100bp single-end sequencing on the Illumina HiSeq 2000 instrument. The sequences obtained were submitted to the Sequence Read Archive at NCBI (PRJNA378655).

## RNA-seq data analysis and functional analysis of differentially expressed genes between control and infected bees

RNA-seq data analysis was performed as reported by Badaoui et al [30]. Briefly, sequence visualization and statistical analyses were performed with FASTQC (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqch). Reads were mapped to the latest version of the honeybee reference genome (Amel_4.5_scaffolds.fa) [22]) using TopHat v2.0.10 with default parameters [31]. Transcript assembly was performed using Cufflinks v2.2.0 [32]. Once the short read sequences were assembled, the output files were merged using the Cuffmerge [31] function and processed using Cuffcompare [32], along with a reference GTF annotation file downloaded from the bee-base database (http://hymenopteragenome.org/beebase/). Differential expression analysis was performed using the Cuffdiff function [32]. Three comparisons of gene expression were made between the control and infected bees at days 5, 10 and 15 P.I. For all the analyses, a feature was considered significant if the false discovery rate (FDR) was less than 0.05. The differentially expressed genes at 5, 10 and 15 days P.I were mapped to pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [33]. The statistical analyses for bee survival and hunger level were performed using the R software "stats" package.

## Proboscis Extension Reflex assay (PER) to assess hunger levels by sucrose uptake behavior

After day 14, eight bees from each box were captured individually with forceps and each bee was placed individually in a glass vial and chilled on ice for 1 min. The bees were fixed to a long plastic drinking straw with tape around the thorax. Testing began 1h17min after the last bees were attached to the straws. The antennae were touched with a droplet of sucrose and if a bee responded by fully extending the proboscis—a Proboscis Extension Response (PER)—was recorded. Each bee was assayed for PER with a concentration series of 0.1%, 0.3%, 1%, 3%, 10% and 30% sucrose solution. Desensitization with water was done between each sucrose concentration. The number of *N. ceranae* spores in each bee was determined under a microscope after the PER test.
## Results

### Honeybees infection with *N. ceranae*: Infection status, survival rates and hunger analysis

Infection with *N. ceranae* was assessed by PCR and qPCR: only three out of five bees sampled from the "inoculated box" at 5 day P.I. were positive for *N. ceranae*, while all 5 bees were positive at days 10 and 15 P.I. *N. ceranae* was not detected in any bees from the "control box". After day 14 P.I., 15 bees were examined for the presence of Nosema spores under the microscope and the number of spores recorded. All the inoculated bees were infected, with a mean of 5.0 x10⁵ spores per bee, and all the control bees were free from any spore (S1 Table).

Survival of infected bees at day 5 P.I. was significantly decreased (p<0.1): 83% compared to 100% for the control bees (S1 Fig). However, the survival rates between the infected and control bees, were not statistically significant at days 10 and 15 P.I. Bees infected with *N. ceranae* had increased levels of hunger compared with controls, assessed by the Proboscis Extension Response and feeding experiments (S2 Fig). An increase in the Proboscis Extension Response was seen especially at intermediate and high sucrose concentrations; 10% of sucrose solution (p-value<0.1) and 30% (p-value<0.1).

## 1. Differentially expressed genes, isoforms and TSSs, and differential promoter usage in honeybees infected with *N. ceranae*

1.1. RNA-seq reads features and statistics. The mean of the reads produced from the eighteen libraries was 13.8 M per library. The mapping rate was, on average, 70% except for one library which had a low mapping rate of 29% and which was discarded from further analyses (S2 Table). The number of expressed unique genes and isoforms in infected and non-infected bees were 15694 and 33182, respectively. The most highly expressed genes, with Fragments Per Kilobase of exon per Million fragments mapped (FPKM) between 1000 and 10000, represented a small fraction of the total genes (0.9%). Genes expressed at a low level, with FPKM between 1 and 10, were considerably more abundant (30–32%). Genes with intermediate expression, FPKM values from 100 to 1000 were the most abundant (33–36%).

The differentially expressed genes, isoforms, TSS and differential promoter usage between control and infected bees were respectively {372;226;274;8} at day 5 P.I., {158;85;116;8} at day 10 P.I. and {67;27;42;0} at day 15 P.I. (S3 Table). Among the differentially expressed genes, 53%, 50% and 30% were also differentially expressed at the isoform level, at days 5, 10 and 15 P.I., respectively (see Venn diagram, S3 Fig). The transcripts with altered forms at more than one time point of infection (days 5, 10 and 15 P.I.) were 15, 7 and 12 for the genes, isoforms and TSS, respectively (see Venn diagram in Fig 1).

1.2. Transcriptional, post-transcriptional regulation and differential promoter usage in bees after infection with *N. ceranae*. Generally, transcriptionally regulated isoforms had different TSSs whilst post-transcriptionally regulated isoforms had the same TSSs (S4 Fig, [24]). The transcriptional start sites and isoforms of the top twenty most differentially expressed genes (S3 Table) at each time point (days 5, 10 and 15 P.I.) were investigated to assess if they were regulated at transcriptional or post-transcriptional levels, or both.

Three groups of genes were defined:

Group1: "Un-spliced and Transcriptionally regulated" genes that had one isoform and one TSS. This group had 9, 6 and 10 genes at days 5, 10 and 15 P.I. respectively (S5, S6 and S7 Figs (group 1)). Most of the genes in this group have a role in the immune response (e.g. GB40713, GB43007, and GB42217)
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Venn diagram illustrating the significantly affected genes (A), isoforms (B) and TSS (C) between infections at days 5, 10 and 15 P.I.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0173438.g001

Group 2: “Spliced and Transcriptionally regulated” genes that had more than one isoform but only one TSS for each isoform. This group had 5, 8 and 3 genes at days 5, 10 and 15 P.I, respectively (S5, S6 and S7 Figs (group 2)). Many of these genes are involved in carbohydrate metabolism and ion transport (e.g. GB51580, GB45152, and GB46516)

Group 3: “Spliced and both Transcriptionally and Post-Transcriptionally regulated” genes with more than one isoform and more than one TSS. This group, contained genes with isoforms that are transcriptionally regulated and other isoforms which are post-transcriptionally regulated and had 6, 6 and 7 genes at days 5, 10 and 15 P.I, respectively (S5, S6 and S7 Figs (group 3)). These genes are involved, among other functions, in cuticle biosynthesis, trans-membrane transport activity and ATP binding (e.g. GB42218, GB45300 and GB47723).

Differential use of promoters between the control and infected bees was investigated as described by Mortazavi *et al* [24]. This approach groups the primary transcripts of a gene based on the promoter used, followed by testing isoform abundance between the control and the infected bees. Eight genes showed different promoter use at day 5 P.I (GB46442, GB41734, GB42277, GB44487, GB46749, GB42364, GB44323 and GB40485), 8 genes at day 10 P.I (GB40787, GB48346, GB47239, GB43964, GB40688, GB42265, GB55747 and GB55915) and none at day 15 P.I (S3 Table).

**1.3. Gene expression characteristics at days 5, 10 and 15 P.I. with N. ceranae.** Principal component analysis (PCA) of the 675 genes which were differentially expressed in at least one pairwise comparison (control vs infection at day 5, control vs infection at day 10, control vs infections at day 15) (Fig 2) showed a clear difference between samples collected at days 5, 10 and 15. The control and infected samples showed the greatest difference at day 5 P.I while this difference was lower at day 10 P.I and no difference between controls and infected bees was found at day 15 P.I. This is reflected in the number of genes that were found to be differentially expressed between the infected bees and controls: 372, 158 and 93 at days 5, 10 and 15 P.I, respectively. The overlap between differentially expressed genes at different time
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Principal component analysis of RNA-seq data. Gene expression changes were investigated at days 5, 10 and 15 P.I. in honeybees infected (INF) with *N. ceranae* or no treatment (CR). The PCA was performed using normalized RNA-Seq data of 675 genes differentially expressed in at least one pairwise comparison: control vs infection at day 5 10 P.I or 15 P.I. Clear differences were seen between samples collected at days 5, 10 and 15 suggesting that ageing of the bees has a larger effect on the pattern of gene expression pattern than infection status.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0173438.g002

points was greater than expected by chance for 5 days P.I versus 10 days P.I (p = 1.77e-52, hypergeometric test), 5 days P.I versus 15 days P.I (p = 6.55e-33) and 10 days P.I versus 15 days P.I (p = 1.92e-33).

At day 5 PI there was a strong down-regulation of Royal jelly genes: *Mrjp1*, *Mrjp2*, *Mrjp3*, *Mrjp4* and *Mrjp6* (S4 Table) and a strong up-regulation of eight serine proteases: *SP22*, *SP40*, *SP44*, *SP17*, *SP18*, *SP35*, *SP36*, *SPH50* as well as the down-regulation of three other serine proteases: *SP34*, *SPH19* (serine protease homolog 19) and *SPH42* (serine protease homolog 42) (S3 and S4 Tables).

**1.4. Pathways affected during honeybee infection by *N. ceranae*.** The most enriched KEGG pathway at all-time points (5, 10 and 15 days P.I) was 'metabolic pathways' (S5A Table, S8 Fig) which includes 'Energy metabolism', 'carbohydrate metabolism', 'amino acids metabolism' and 'lipid metabolism'. The number of genes differentially expressed between infected and control bees in the metabolic pathways decreases considerably from day 5 P.I (29 genes) to day 10 P.I (17 genes) and day 15 P.I (9 genes). The changes in the expression of 29 genes in this pathway were followed through the time-course of the challenge (Fig 3). Except the chitinase 5 (*Cht5*) gene expression of which increased from day 5 P.I to day 10 P.I and to day 15 P.I, expression of all the other genes decreased from day 5 P.I to days 10/15 P.I. Some genes reached their lowest level of expression at day 10 P.I (GB44841, GB45654, GB50218, GB51814, GB42963, GB49240, GB55706, GB42434), while others had the lowest level at day 15 P.I (GB52756, GB51751, GB52724, GB42425, GB52923, GB48905, GB44367, GB40022, GB55705) (Fig 3). The number of differentially expressed genes in the lysosome pathway, which is within the "metabolic pathways", was also enriched at day 5 P.I (GB53306, GB43825, GB51751, GB54097, GB55242, GB42616) and day 10 P.I (GB53306, GB43825, GB51751, GB53579) although not at day 15 P.I. (S5A Table).
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Expression profiles of the genes involved in the 'metabolic pathways'. The X axis shows the time points following infection (days 5, 10 and 15 P.I). On the Y axis, genes expression fold changes.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0173438.g003

The KEGG annotation database contains little information on the bee signaling pathways. Therefore the differentially expressed genes were used to manually infer other pathways potentially modified during infection by *N. ceranae*. This identified the 'Immunity pathway', 'Trehalose transporter pathway', 'Cuticular pathway', 'Cytochrome pathway', 'Serine protease pathway' and 'Neurogenesis pathway' as being affected by the infection (S5B Table).

## Discussion

In this work we report differentially expressed genes, isoforms, TSS and promoter usage associated with infection of the bees with *N. ceranae*.

## Regulation complexity of differentially expressed genes in bees infected with *N. ceranae*

The most differentially expressed genes between uninfected and infected bees were classified into three groups according to their transcriptional/post-transcriptional regulation. The first two groups included genes involved in the innate immune response. These genes did not show any evidence of different isoforms and/or post-transcriptional regulation. Examples of these groups include: *GB48966*, a homeotic protein caudal, which operates post-embryonically in the intestine where it regulates antimicrobial peptide levels to maintain the normal gut flora [34]; *GB51833*, Sodium dependent nutrient amino acid trasporter, which contributes to the synthesis of catecholamines for sclerotization pathways involved incuticle colour, wound curing, and immune responses [35].
Post transcriptional regulation is responsible for expanding protein functional diversity and differs dependent on environmental conditions and developmental stage to expand phenotypic plasticity on the molecular level [36,37]. However, the selection of the incorrect splices sites may cause gene dysfunction [38]. The genes in these first two groups may be under strong selection pressure to prevent the emergence of new isoforms or alternative transcriptional regulation to avoid loss of function.

*N. ceranae* is amitochondriate, and so has high dependency on host ATP. The energy demands of the parasite increase the energy requirement for the bee, which is seen in the increased hunger of infected bees revealed by the PER test. GB41033, which is involved in the metabolism of the sugar, was highly expressed in infected bees, and is consistent with increased sugar metabolism associated with the higher energy requirement [11]. *GB42217*, *acyl-CoA desaturase 1-like*, has a role in oxidation-reduction process following infection [39].

The third group includes genes mainly regulate biochemical and developmental processes, and includes receptors such as serpentine (GB45151) that regulate the expression of multiple classes of developmental genes [39, 40], SLC6 (GB51834) that transport a whole suite of molecules across cell membranes [41] and hedgehog (GB45300) which is required for normal regulation of cell proliferation and differentiation during embryonic stages [40].

This third group genes includes immune function genes and shows abundant examples of alternative splicing. There are many examples where the immune system has been shown to use diverse mechanisms of gene regulation to increase the versatility of response, including alternative splicing, different TSS and promoter usage [42,43], for example, in vertebrates alternative splicing of CD44, a protein involved in T cell homing, is crucial for T cell function [43]. The transcription factors p63, p73 and p53, implicated in cell response to stress and development, and encodes multiple proteins as a result of post transcriptional regulation [42]. Diverse post-transcriptional regulation mechanisms seen for genes involved in bee immune response, cell signals and development may increase the immune repertoire to respond to infection.

Alternative promoter usage was found for 8 genes at days 5 and 10P.I in comparison with uninfected controls, but none were seen at day 15 P.I. This suggests that alternative promoter usage is important for bee-*N. ceranae* interaction as seen in previous studies where differential promoter usage has been associated with the diversity of gene expression pattern to deal with different pathogen challenges [44,45]. That no alternative promoter usage was found at day 15 P.I and that a small number of differentially expressed genes is seen at this time point compared with time points 5 and 10 P.I. may suggest that the host and parasite are becoming adjusted to each other. It has been reported that genes with alternative promoters use are also more likely to be differentially expressed [46]. Genes related to human cancer which have altered expression have on average 2 promoters compared with an average of 1.5 promoters for the other human genes [47].

The overrepresentation of the splice variants, 30%-53% of genes between 5 and 10 days PI respectively of differentially expressed genes, seen following infection, suggests that extending protein diversity is important for bees to deal with the *N. ceranea* infection. Therefore, start site and post-transcriptional regulation is potentially important for tuning the response of the bee to *N. ceranea* infection, possibly by controlling the strength and the duration of the response [48]. This is in accordance with findings reported for Drosophila melanogaster [49] and mammals [50] where splice variants play an important role in dealing with infection. Interestingly, alternative TSSs have been associated with stable differences in behavior in honeybees [51].

Genes involved in the “metabolic pathways” (S5 Table) have higher expression at day 5 P.I than days 10 and 15 P.I. This suggests an increased energy demand placed on the bee by *N. ceranae* during early stages of the infection which is in agreement with the PER test responses.
Pathways enriched in bees infected with *N. ceranae*

The global pathway 'metabolic pathways' showed a general increase in expression at day 5 P.I. Increased energy metabolism via carbohydrate catabolism is a key feature of bees infected with *N. ceranae*. Increased energy need and hence hunger was seen in the present study from the PER tests. *Cht5*, *GMCOX3*, *LOC55124* and *LOC408474* were the most highly expressed genes in this pathway at day 5 P.I with a log fold change greater than 1.5. *Cht5* and *GMCOX3* are involved in chitin metabolism and immune response respectively, and have previously been found to be modulated during honeybee infection with *N. cerania* [25]. While *GMCOX3* expression decreased after day 5 P.I, expression of *Cht5* continued to increase to reach log fold changes of 4.28 and 4.52 at days 5 and 10 P.I, respectively. *LOC55124* and *LOC408474* are involved in protein and lipid metabolism, respectively.

The lysosome pathway also plays an important role in bees infected with *N. ceranae* which is shown by the genes GB53306, GB43825, GB51751, GB54097, GB55242, and GB42616in this pathway having increased levels of expression in infected bees. Lysosome has been implicated in antifungal and antiviral activity [52]. An increased level of lysosome activity maybe a defense strategy used by the bee to digest the *N. ceranae* mycelia. Furthermore, the pathways: 'Immunity', 'Trehalose transporter', 'Cuticular', 'Cytochrome', 'Serine protease' and 'Neurogenesis' were modulated in infected bees compared with controls (S5B Table). These pathways have also been found to be regulated by bees in response to *N. ceranae*, *N. Microsporidia*, *E.coli* and Paenibacillus larva infections [25, 11, 53].

Infection with *N. ceranae* reduces the bees energy reserves

Food quantity and quality is a key factor regulating the development of female larvae, either into a queen or a worker. Specifically, the Major Royal Jelly Proteins (MRJPs) are consumed to some extent by all larvae, but mainly by the larvae destine to become the adult queen [35]. The Royal Jelly genes: *Mrjp1*, *Mrjp2*, *Mrjp3*, *Mrjp4* and *Mrjp6* were down regulated at day 5 P.I. The protein hexamerin (*HEX*) is accumulated by the bees to enhance their growth [54]. Expression of this gene was not changed at day 5 P.I but was down-regulated at days 10 and 15 P.I. Down-regulation of *HEX110*, *HEX70b* and *HEX70c* has previously been reported in bees infected with *Paenibacillus larvae* and *chalkbrood fungus* [55, 56]. Neuropeptide Y (NPY), has been implicated in regulating food intake [57]. The receptor for this neuropeptide was up-regulated following *N. ceranae* infection. Insulin growth factor binding protein was also up-regulated (S3 Table). Insulin signaling regulates nutritional balance [57]. Importantly, among the top pathways differentially expressed (up-regulated) at day 5 P.I, were 'Energy metabolism', 'lipid metabolism' and 'carbohydrate metabolism'. Therefore *N. ceranae* infection seems to affect metabolism and nutritional status of the infected bees. The increased transcription of four genes mapping to the trehalose transporter (*LOC413576*, *LOC726505*, *LOC724874* and *LOC413575*), which constitute the main major carbohydrate energy storage molecule in insects, supports this observation.

Many studies have reported changes in nutrition [14,58] and increased feeding rate [59,60] in parasitized insects. In particular, studies of bees infected with *N. ceranae* have reported a higher sugar demand [11,14,60–64]. *N. ceranae* requires ATP to replicate within the honey bee mid-gut tissue [11,65] but being amitochondriate, has limited metabolic capacity and is unable to produce ATP. As a consequence it is completely dependent on host ATP [66,67]. *N. ceranae* therefore imposes an energetic demand on its host, which is consistent with the changes expression of metabolism related genes the increased hunger responses and lower survival seen in the present study.
*N. ceranae* infection alters amino acid metabolism and suppresses the bee immune response

The amino acids metabolic pathway was modulated by *N. ceranae* infection. In a previous study, worker bees infected with *N. ceranae* were found to have altered amino acids titers in their haemolymph [68], which may affect the immune function [69–71]. Disturbance of protein metabolism may be a strategy adopted by *N. ceranae* to suppress the bee's immune response [13] and promote its own survival. Bees infected with Varroa mite, show changes the expression of genes related to protein metabolism, including the down regulation of the Vg [72–74]. In the present study the expression of genes coding for the Vg, MRJP protein were found to be down-regulated.

Following infection by *N. ceranae*, a large number of genes coding for serine proteases were modulated (*SP17*, *SP18*, *SP22*, *SP34*, *SP35*, *SP36*, *SP40*, *SP44*, *SPH19*, *SPH42*,*SPH50*) at days 5, 10 and 15 P.I. The serine proteases are implicated in the immune response for Drosophila where 94 out of 201 serine protease genes are involved in diverse immune protease cascades [53].

Except SP36 and SPH42, all these genes showed alternative splicing and differential TSS (S3 Table) which might increase their expression plasticity and function. Activation of these molecules in the arthropod hemolymph is a central component of several immune responses, including the activation of antimicrobial peptide synthesis, and modulation of hemocyte function [75, 76].

Therefore, in addition to placing additional energy demands on the host *N. ceranae* also seems to affect the immune response of infected bees.

Honeybees have many immune pathways and defense mechanisms consisting of both cellular and humoral immune response [16] including antimicrobial peptide (AMPs) synthesis [77]. Honeybee humoral immunity involves four main antimicrobial peptides: abaecin, apidaecin, defensin and hymenoptaecin [20]. In the present study, apidaecin (*Apid1*, *Apid73*), defencin (*Def1*) and hymenoptaecin (*GB51223*) were strongly down-regulated at day 5 P.I compared with uninfected controls but not at days 10 and 15 P.I, except for the hymenoptaecin gene that was up-regulated at day 15 P.I. All these genes have at least two differentially expressed isoforms and TSS [78]. Different sequences of apidaecin arising from post translational regulation have been reported [79, 17], and Klaudiny et al [78] have reported a novel defensin isoform in honeybee. The finding reported here suggests that specific splice variation and TSS may occur in response to pathogens.

The egg yolk protein vitellogenin, which has been implicated in regulating the number of haemocytes [80] was significantly decreased following *N. ceranae* infection at day 5 P.I. Immune suppression has been reported following infection with by *N. ceranae* [13] and the Varroa mite [81].

Two genes coding for the major components of pathogen recognition (PGRPs) were modulated in this study: peptidoglycan recognition protein S1 (*PGRP-S1*) was down-regulated and peptidoglycan recognition protein S3 (*PGRP-S3*) was up-regulated. Honeybees have four PGRPs, PGRP-S1, PGRP-S2, PGRP-S3, and PGRP-LC [79], either act as recognition proteins, or can degrade bacterial cell wall through amidase activity [82, 83]. Honeybees have much lower diversity of PGRP than Drosophila and Anopheles which have 13 and seven genes, respectively. This may limit the diversity of pathogens that bees can recognize compared with drosophila [79]. The present study identified that all the honeybee PGRPs were subjected to alternative splicing and TSS differentiation, which may be an evolutionary strategy to diversify the protein repertoire and overcome the limited number of PGRPs.
Cuticular, olfactory and neuronal modifications during infection of bees
with *N. ceranae*

Cuticle coatings are important to protect insects from pathogens and cuticle structure may be altered during infection [84]. In this work, seven cuticle genes (LOC724464, CPR19, LOC725547, Cpap3-a, CPR24, CPR5 and Cpap3-dC) were down-regulated at day 5 P.I compared with controls and one additional gene (CPR2) was down-regulated at day 10 P.I. All these genes showed alternative splicing and differential TSS. The down-regulation of these cuticle genes maybe a strategy of the pathogen to improve its survival, by facilitating transfer among members of infected hive. Expression of the cuticle genes is associated with hygienic behavior [85, 86] which may also affect parasite transfer and survival.

Infection by *N. ceranae* may also change the behavior of bees by altering olfactory acuity. Three genes coding for (*Obp14*, *Obp18* and *Obp3*), were down-regulated at day 5 P.I. compared with controls, two (*obp3*, *obp17*) at day 10 P.I and two (*obp14*, *obp3*) at day 15 P.I. The odor binding proteins are not restricted to olfaction and they have other functions, including the adaptation to different environments [87]. Changes in expression of odor binding proteins following *N. ceranae* may change the behavioral or physiological responses of bees.

Transcripts involved in neurogenesis such as Slit genes, LOC410555, LOC724772, leucine-rich repeat neuronal protein 3 (LOC724187), neuropeptide and neurotransmitter transporter 8, were modulated at day 5 P.I but not at days 10 and 15 P.I. Slitis part of a large network of genes involved in tissue regeneration and has been shown to be down regulated following infection by Nosema [11]. Changes in expression of genes involved in neuronal development have also been reported following infestation with the Varroa mite [88]. Therefore bee parasites may affect behavior either to enhance their survival or to promote the spread within the colony.

Conclusions

*N. ceranae* infection results in the regulation of gene function, both at the level of expression and through changes in promoter usage and post translational variation. The effects of *N. cera-nae* infection are seen on expression of genes involved in three key processes:1) energetic stress which is reflected in elevated hunger levels, 2) the modulation of the metabolic pathways, par-ticularly immune function and 3) changes in behavior mediated through altered olfactory or hygiene related gene expression. These observed changes could be associated with the host responding to clear the infection, and with the infecting *N.Ceranea* altering host gene expres-sion to promote its survival.

**Supporting information**

S1 Table. Counts of *N. ceranae* in control and infected bees with *N.Ceranea*. 100% of the inoculated bees used for the bioassay were infected by Nosema while no control bees were infected.
(PDF)

S2 Table. (i) RNA-seq reads numbers generated per sample, (ii) RNA-Seq reads number and percentage mapped to the honeybee genome and (iii) the number and percentage of RNA-seq reads that mapped to more than one site in the genome.
(PDF)

S3 Table. Differentially expressed genes, isoforms, TSS and differential promoter usage during infection with *N.Ceranea* at days 5, 10 and 15 P.I.
(XLS)
**S4 Table.** Genes differentially expressed (|FC|> = 1.5) between infected and control bees at days 5, 10 and 15 P.I.
(PDF)

**S5 Table.** A. mapping of the differentially expressed genes at days 5, 10 and 15 P.I. to pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The top pathways in which at least 4 genes were involved were selected. B. Manually inferred pathways using the differentially expressed genes at days 5, 10 and 15 P.I.
(ODS)

**S1 Fig.** Survival rates of control and infected bees with *N. ceranae* up to 14 days Post Infec-
tion (P.I.).
(PDF)

**S2 Fig.** Percentage of bees triggering a Proboscis Extension Response at different concen-
tration sucrose solution (30%,10%, 3%, 1%, 0.3% and 0.1%).
(PDF)

**S3 Fig.** Venn diagram illustrating the significantly affected (A) genes (B), isoforms and (C) TSSs found in bees infected with *N. ceranae* at days 5, 10 and 15 P.I.
(PDF)

**S4 Fig.** Cufflinks approach for estimating transcriptional and post-transcriptional regula-
tory effects on overall transcript output. Let's consider a transcript with two isoforms
(A, B). When the abundance of isoforms A, B and C are grouped by TSS, the changes in the
relative abundance of the TSS groups indicate transcriptional regulation (A+B vs. C). Post-
transcriptional effects are observed as changes in the levels of the isoforms in a single TSS
group (A vs. B) (Adapted from Trapnell et al. 2012).
(PDF)

**S5 Fig.** Transcriptional/Post-transcriptional regulation of the top twenty significant genes following the infection of bees with *N. ceranae* at day 5 P.I. (Group 1) Un-spliced and tran-
scriptionally regulated genes,(Group 2) spliced and post-transcriptionally regulated genes and
(Group 3) spliced and both transcriptionally and post-transcriptionally regulated genes. For
each transcript, the "XLOC", "TSS" and "TCONS" suffixes correspond to the genes, TSSs and
isoforms, respectively. Differentially expressed isoforms with different TSSs are transcription-
ally regulated, while isoforms with the same TSS are regulated at the post-transcriptional level.
(PDF)

**S6 Fig.** Transcriptional/Post-transcriptional regulation of the top twenty significant genes following the infection of bees with *N. ceranae* at day 10 P.I. (Group 1) Un-spliced and tran-
scriptionally regulated genes,(Group 2) spliced and post-transcriptionally regulated genes and
(Group 3) spliced and both transcriptionally and post-transcriptionally regulated genes. For
each transcript,the "XLOC", "TSS" and "TCONS" suffixes correspond to the genes, TSSs and
isoforms, respectively. Differentially expressed isoforms with different TSSs are transcription-
ally regulated, while isoforms with the same TSS are regulated at the post-transcriptional level.
(PDF)

**S7 Fig.** Transcriptional/Post-transcriptional regulation of the top twenty significant genes following the infection of bees with *N. ceranae* at day 15 P.I. (Group 1) Un-spliced and tran-
scriptionally regulated genes,(Group 2) spliced and post-transcriptionally regulated genes and
(Group 3) spliced and both transcriptionally and post-transcriptionally regulated genes. For
each transcript, the "XLOC", "TSS" and "TCONS" suffixes correspond to the genes, TSSs and
isoforms, respectively. Differentially expressed isoforms with different TSSs are transcription-
ally regulated, while isoforms with the same TSS are regulated at the post-transcriptional level.
(PDF)

**S8 Fig. The 'metabolic pathways' enriched following the infection of bees with *N.Ceranea*.**
This global pathway includes mainly 'energy metabolism', 'carbohydrate metabolism', 'amino acids metabolism' and 'lipid metabolism' related genes.
(TIF)

**Acknowledgments**

This project was supported by the Eureka-Eurostars Project E!5069 NewBone. The authors
would like to thank Maria Cecere for the technical help. The authors declare no competing
financial interests.

**Author Contributions**

**Conceptualization:** BB SB EG.

**Data curation:** BB AA MC.

**Formal analysis:** BB SB AF FP CG AC AG GC GF FM.

**Funding acquisition:** EG SB.

**Investigation:** BB SB AF FP CG AC AG GC GF FM.

**Methodology:** BB SB AF FP CG AC AG GC GF FM EG.

**Project administration:** EG SB BB.

**Resources:** SB AF FP CG AC AG GC GF FM.

**Software:** BB AA MC.

**Supervision:** SB.

**Validation:** BB SB AF FP CG AC AG GC GF FM.

**Visualization:** BB JLW.

**Writing – original draft:** BB SB.

**Writing – review & editing:** BB JLW.

**References**

1. Higes M, Meana A, Bartolome C, Botias C, Martin-Hernandez R. Nosema ceranae (Microsporidia), a controversial 21st century honey bee pathogen. Environ Microbiol Rep. 2013; 5(1):17–29. https://doi.org/10.1111/1758-2229.12024 PMID: 23757127

2. Paxton RJ, Klee J, Korpela S, Fries I. Nosema ceranae has infected Apis mellifera in Europe since at least 1998 and may be more virulent than Nosema apis. Apidologie. 2007; 38(6):558–65.

3. Higes M, Martín-Hernández R, Garrido-Bailà E, Botias C, Meana An. The presence of Nosema ceranae (Microsporidia) in North African honey bees (Apis mellifera intermissa). Journal of Apicultural Research. 2009; 48(3):217–9.

4. Martín-Hernández R, Meana A, Prieto L, Salvador AM, Garrido-Bailòn E, Higes M. Outcome of colonization of Apis mellifera by Nosema ceranae. Applied and Environmental Microbiology. 2007; 73 (20):6331–8. https://doi.org/10.1128/AEM.00270-07 PMID: 17675417
5. Higes M, Martin-Hernández R, Botías C, Bailàn EG, González-Porto AV, Barrios L, et al. How natural infection by Nosema ceranae causes honeybee colony collapse. Environmental microbiology. 2008; 10 (10):2659–69. https://doi.org/10.1111/j.1462-2920.2008.01687.x PMID: 18647336

6. Topolska Gy, Gajda A, Hartwig A. Polish honey bee colony-loss during the winter of 2007/2008. Journal of Apicultural Science. 2008; 2(52) 95–103.

7. Bacandritsos N, Granato A, Budge G, Papanastasiou I, Roinioti E, Caldon M, et al. Sudden deaths and colony population decline in Greek honey bee colonies. J Invertebr Pathol. 2010; 105(3):335–40. https://doi.org/10.1016/j.jip.2010.08.004 PMID: 20804765

8. Borneck R, Viry A, Martin-Hernández R, Higes M. Honey bee colony losses in the Jura Region, France and related pathogens. Journal of Apicultural Research. 2010; 49(4):334–6.

9. Currie RW, Pernal SF, Guzmàn-Novoa E. Honey bee colony losses in Canada. Journal of Apicultural Research. 2010; 49(1):104–6.

10. Van der Zee Re. Colony losses in the Netherlands. Journal of Apicultural Research. 2010; 49(1):121–3.

11. Dussaubat C, Brunet JL, Higes M, Colbourne JK, Lopez J, Choi JH, et al. Gut pathology and responses to the microsporidium Nosema ceranae in the honey bee Apis mellifera. PLoS One. 2012; 7(5):e37017. https://doi.org/10.1371/journal.pone.0037017 PMID: 22623972

12. Garcia-Palencia P, Martin-Hernandez R, Gonzalez-Porto A-V, Marin P, Meana A, Higes M. Natural infection by Nosema ceranae causes similar lesions as in experimentally infected caged-worker honey bees (Apis mellifera). Journal of Apicultural Research. 2010; 49(3):278–83.

13. Antunez K, Martin-Hernández R, Prieto L, Meana A, Zunino P, Higes M. Immune suppression in the honey bee (Apis mellifera) following infection by Nosema ceranae (Microsporidia). Environ Microbiol. 2009; 11(9):2284–90. https://doi.org/10.1111/j.1462-2920.2009.01953.x PMID: 19737304

14. Mayack C, Naug D. Energetic stress in the honeybee Apis mellifera from Nosema ceranae infection. J Invertebr Pathol. 2009; 100(3):185–8. https://doi.org/10.1016/j.jip.2008.12.001 PMID: 19135448

15. Aliferis KA, Copley T, Jabaji S. Gas chromatography-mass spectrometry metabolite profiling of worker honey bee (Apis mellifera L.) hemolymph for the study of Nosema ceranae infection. J Insect Physiol. 2012; 58(10):1349–59. https://doi.org/10.1016/j.jinsphys.2012.07.010 PMID: 22841888

16. Osta MA, Christophides GK, Vlachou D, Kafatos FC. Innate immunity in the malaria vector Anopheles gambiae: comparative and functional genomics. J Exp Biol. 2004; 207(15):2551–63.

17. Casteels P, Ampe C, Jacobs F, Vaeck M, Tempst P. Apidaecins: antibacterial peptides from honeybees. EMBO J. 1989; 8(8):2387–91. PMID: 2676519

18. Casteels P, Ampe C, Riviere L, Van Damme J, Elicone C, Fleming M, et al. Isolation and characterization of abaecin, a major antibacterial response peptide in the honeybee (Apis mellifera). Eur J Biochem. 1990; 187(2):381–6. PMID: 2298215

19. Casteels P, Ampe C, Jacobs F, Tempst P. Functional and chemical characterization of Hymenoptaecin, an antibacterial polypeptide that is infection-inducible in the honeybee (Apis mellifera). Journal of Biological Chemistry. 1993; 268(10):7044–54. PMID: 8463238

20. Casteels-Josson K, Zhang W, Capaci T, Casteels P, Tempst P. Acute transcriptional response of the honeybee peptide-antibiotics gene repertoire and required post-translational conversion of the precursor structures. J Biol Chem. 1994; 269(46):28569–75. PMID: 7961803

21. Consortium HBGS. Insights into social insects from the genome of the honeybee Apis mellifera. Nature. 2006; 443(7114):931–49. https://doi.org/10.1038/nature05260 PMID: 17073008

22. Elsik CG, Worley KC, Bennett AK, Beye M, Camara F, Childers CP, et al. Finding the missing honey bee genes: lessons learned from a genome upgrade. BMC genomics. 2014; 15(1):86.

23. Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, Gerstein M, et al. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science. 2008; 320(5881):1344–9. https://doi.org/10.1126/science.1158441 PMID: 18451266

24. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptiones by RNA-Seq. Nat Methods. 2008; 5(7):621–8. https://doi.org/10.1038/nmeth.1226 PMID: 18516045

25. Aufauvre J, Misme-Aucouturier B, Vigues B, Texier C, Delbac F, Blot N. Transcriptome analyses of the honeybee response to Nosema ceranae and insecticides. PLoS One. 2014; 9(3):e91686. https://doi.org/10.1371/journal.pone.0091686 PMID: 24646894

26. McDonnell CM, Alaux C, Parrinello H, Desvignes JP, Crauser D, Durbesson E, et al. Ecto- and endoparasite induce similar chemical and brain neurogenomic responses in the honey bee (Apis mellifera). BMC Ecol. 2013; 13:25. https://doi.org/10.1186/1472-6785-13-25 PMID: 23866001

27. Ward KN, Coleman JL, Clinnin K, Fahrbach S, Rueppell O. Age, caste, and behavior determine the replicative activity of intestinal stem cells in honeybees (Apis mellifera L.). Experimental Gerontology 2008; 43:530–537 https://doi.org/10.1016/j.exger.2008.03.012 PMID: 18479865
28. Kux K, Pitsouli C Tissue communication in regenerative inflammatory signaling: lessons from the fly gut. Front. Cell Infect. Microbiol. 2014;472-78.

29. Chen Y, Evans JD, Zhou L, Boncristiani H, Kimura K, Xiao T, et al. Asymmetrical coexistence of Nosema ceranae and Nosema apis in honey bees. J Invertebr Pathol. 2009; 101(3):204-9. https://doi.org/10.1016/j.jip.2009.05.012 PMID: 19467238

30. Badaoui B, Rutigliano T, Anselmo A, Vanhee M, Nauwynck H, Giuffra E, et al. RNA-sequence analysis of primary alveolar macrophages after in vitro infection with porcine reproductive and respiratory syndrome virus strains of differing virulence. PloS one. 2014; 9(3).

31. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009; 10(3):R25. https://doi.org/10.1186/gb-2009-10-3-r25 PMID: 19261174

32. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation.2010; Nature biotechnology. 28(5):511-5. https://doi.org/10.1038/nbt.1621 PMID: 20436464

33. Kanehisa M et al. Data, information, knowledge and principle: back to metabolism in KEGG. Nucleic Acids Res. 42. 2014; D199–D205 https://doi.org/10.1093/nar/gkt1076 PMID: 24214961

34. Hwang M.S. et al. Nagy L.M. et al. Schulz C. et al. Moreno E. et al. Shinmyo Y. et al. The caudal homeodomain protein activates Drosophila E2F gene expression. Nucleic Acids Res. 2002; 30:5029-5035. PMID: 12466526

35. Snodgrass RE. Anatomy of the honey bee: Cornell University Press; 1984.

36. Mastrangelo AM, Marone D, Laidò G, De Leonardis AM, and De Vita P. Alternative splicing: enhancing ability to cope with stress via transcriptome plasticity. Plant Sci. 2012; 185-186: 40-49.

37. Chen J, Nolte V, Schlötterer C, Temperature-Related Reaction Norms of Gene Expression: Regulatory Architecture and Functional Implications. Mol. Biol. Evol. 2015; 32: 2393-2402. https://doi.org/10.1093/molbev/msv120 PMID: 25976350

38. Krawczak M, Thomas NS, Hundrieser B, et al. Single base-pair substitutions in exon—intron junctions of human genes: nature, distribution, and consequences for mRNA splicing. Hum Mutat. 2007; 28:150–8. https://doi.org/10.1002/humu.20400 PMID: 17001642

39. Weinstock et al. Insights into social insects from the genome of the honeybee Apis mellifera. Nature. 2006; 443, 931–949. https://doi.org/10.1038/nature05260 PMID: 17073008

40. Alcedo J, Ayzenzon M, Von Ohlen T, Noll M, Hooper J E. The *Drosophila smoothened* gene encodes a seven-pass membrane protein, a putative receptor for the hedgehog signal. Cell.1996; 86:221-232. PMID: 8706127

41. Saier MH, Yen MR, Noto K, Tamang DG, Elkan C. The Transporter Classification Database: recent advances. Nucleic Acids Res. 2009; 37:D274–D278. https://doi.org/10.1093/nar/gkn862 PMID: 19022853

42. Murray-Zmijewski F AU, -Lane, D P AU, Bourdon, J-C TI. p53//p63//p73 isoforms: an orchestra of isoforms to harmonise cell differentiation and response to stress JA—Cell Death Differ PY. 2006/04/07/ online VL—13 IS—6 SP—962 EP—972 SN—1350–9047.

43. Bas J, Baaten G, Cheng-Rui Li and Bradley Linda M. Multifaceted regulation of T cells by CD44. Communicative & Integrative Biology.2010; 3:6, 508–512.

44. Carninci P, Sandelin A, Lenhard B, Katayama S, Shimokawa K, Ponjavic J, Semple CAM, Taylor MS, Engstrom PG, Frith MC, et al. Genome-wide analysis of mammalian promoter architecture and evolution. Nat Genet.2006; 38: 626–635. https://doi.org/10.1038/ng1789 PMID: 16645617

45. Suzuki H, Forrest ARR, vanNimwegen E, Daub CO, Balwierz PJ, Irvine KM, Lassmann T, Ravasi T, Hasegawa Y, de Hoon MJL, et al. The transcriptional network that controls growth arrest and differentiation in a human myeloid leukemia cell line. Nat Genet. 2009; 41: 553–562. https://doi.org/10.1038/ng.375 PMID: 19377474

46. Liu S. Increasing Alternative Promoter Repertories Is Positively Associated with Differential Expression and Disease Susceptibility. 2010; PLOS ONE 5(3): e9482. https://doi.org/10.1371/journal.pone.0009482 PMID: 20208995

47. Davuluri RV, Suzuki Y, Sugano S, Plass C, Huang TH; The functional consequences of alternative promoter use in mammalian genomes. Trends Genet. 2008; 24: 167-177. https://doi.org/10.1016/j.tig.2008.01.008 PMID: 18329129

48. Modrek B, Lee C. A genomic view of alternative splicing. Nat Genet. 2002; 30: 13–19 https://doi.org/10.1038/ng0102-13 PMID: 11753382

49. Brown J. B. et al. Diversity and dynamics of the Drosophila transcriptome. Nature. 2014; 512, 393-399. https://doi.org/10.1038/nature12962 PMID: 24670639
50. Kawaji H. et al. Dynamic usage of transcription start sites within core promoters. Genome Biol. 7, R118. https://doi.org/10.1186/gb-2006-7-12-r118 PMID: 17156492

51. Khamis AM, Hamilton AR, Medvedeva YA, et al. Insights into the Transcriptional Architecture of Behavioral Plasticity in the Honey Bee *Apis mellifera*. *Scientificreport*.2015; 5:11136.

52. Samaranayake YH, Samaranayake LP, Pow EHN, Beena VT, Yeung KWS. Antifungal Effects of Lysozyme and Lactoferrin against Genetically Similar, Sequential Candida albicans Isolates from a Human Immunodeficiency Virus-Infected Southern Chinese Cohort J Clinic Microbiol.2001; 39:3296-3302

53. Zou Z, Lopez DL, Kanost MR, Evans JD, Jiang H. Comparative analysis of serine protease-related genes in the honey bee genome: possible involvement in embryonic development and innate immunity. Insect Mol. Biol.2006; 15, 603–614. https://doi.org/10.1111/j.1365-2583.2006.00684.x PMID: 17069636

54. Queenie WTC and Leonard JF. Changes in protein expression during honey bee larval development. Genome Biol. 2008; 9(10):R156. https://doi.org/10.1186/gb-2008-9-10-r156 PMID: 18959778

55. Chan QW, Melathopoulos AP, Pernal SF, Foster LJ. The innate immune and systemic response in honey bees to a bacterial pathogen, Paenibacillus larvae. BMC Genomics. 2009; 10:387. https://doi.org/10.1186/1471-2164-10-387 PMID: 19695106

56. Aronstein KA, Murray KD, Saldivar E. Transcriptional responses in honey bee larvae infected with chalkbrood fungus. Biomed Central (BMC) Genomics. 2010. 11:391.

57. Ament SA, Wang Y, Robinson GE. Nutritional regulation of division of labor in honey bees: toward a systems biology perspective. Wiley Interdiscip Rev Syst Biol Med. 2010; 2(5):566-76. https://doi.org/10.1002/wsbm.73 PMID: 20836048

58. Thompson SN, Redak RA. Parasitism of an insect Manduca sexta L. alters feeding behaviour and nutrient utilization to influence developmental success of a parasitoid. J Comp Physiol B. 2008; 178(4):515-27. https://doi.org/10.1007/s00360-007-0244-6 PMID: 18196249

59. Grimstad PR, Ross QE, Craig GB Jr. Aedes triseriatus (Diptera: Culicidae) and La Crosse virus. II. Modification of mosquito feeding behavior by virus infection. J Med Entomol. 1980; 17(1):1-7. PMID: 7365753

60. Holt H, Aronstein K, Grozinger C. Chronic parasitization by Nosema microsporidia causes global expression changes in core nutritional, metabolic and behavioral pathways in honey bee workers (Apis mellifera). BMC Genomics. 2013; 14, 799 https://doi.org/10.1186/1471-2164-14-799 PMID: 24245482

61. Dussaubat C, Brunet J-L, Higes M, Colbourne JK, Lopez J, et al. Gut pathology and responses to the microsporidium Nosema ceranae in the honey bee Apis mellifera. PLoS ONE. 2012. 7: e37017. https://doi.org/10.1371/journal.pone.0037017 PMID: 22623972

62. Naug D, Gibbs A. Behavioral changes mediated by hunger in honeybees infected with Nosema ceranae. Apidologie. 2009; 40(6):595-9.

63. Martin-Hernandez R, Botías C, Barrios L, Martinez-Salvador A, Meana A, Mayack C, et al. Comparison of the energetic stress associated with experimental Nosema ceranae and Nosema apis infection of honeybees (Apis mellifera). Parasitol Res. 2011; 109(3):605-12. https://doi.org/10.1007/s00436-011-2292-9 PMID: 21360094

64. Goncalves RL, Oliveira JH, Oliveira GA, Andersen JF, Oliveira MF, Oliveira PL, et al. Mitochondrial reactive oxygen species modulate mosquito susceptibility to Plasmodium infection. PLoS One. 2012; 7(7):e41083. https://doi.org/10.1371/journal.pone.0041083 PMID: 22815925

65. Malone LA, Gatehouse HS. Effects of Nosema apis Infection on Honey Bee (Apis mellifera) Digestive Proteolytic Enzyme Activity. Journal of Invertebrate Pathology. 1998; 71(2):169-74.

66. Keeling PJ, Fast NM. Microsporidia: biology and evolution of highly reduced intracellular parasites. Annu Rev Microbiol. 2002; 56:93-116. https://doi.org/10.1146/annurev.micro.56.012302.160854 PMID: 12142484

67. Cornman RS, Chen YP, Schatz MC, Street C, Zhao Y, Desany B, et al. Genomic analyses of the microsporidian Nosema ceranae, an emergent pathogen of honey bees. PLoS Pathog. 2009; 5(6):e1000466. https://doi.org/10.1371/journal.ppat.1000466 PMID: 19503607

68. Wang D-I, Moeller FE. Comparison of the free amino acid composition in the hemolymph of healthy and Nosema-infected female honey bees. Journal of invertebrate pathology. 1970; 15(2):202-6.

69. Huang Z. Pollen nutrition affects honey bee stress resistance. Terrestrial Arthropod Reviews. 2012; 5(2):175-89.

70. Lee KP, Cory JS, Wilson K, Raubenheimer D, Simpson SJ. Flexible diet choice offsets protein costs of pathogen resistance in a caterpillar. Proc Biol Sci. 2006; 273(1588):823-9. https://doi.org/10.1098/rspb.2005.3385 PMID: 16618675

71. Ponton F, Wilson K, Holmes AJ, Cotter SC, Raubenheimer D, Simpson SJ. Integrating nutrition and immunology: a new frontier. J Insect Physiol. 2013; 59(2):130-7. https://doi.org/10.1016/j.jinsphys. 2012.10.011 PMID: 23159523
72. Alaux C, Dantec C, Parrinello H, Le Conte Y. Nutrigenomics in honey bees: digital gene expression analysis of pollen's nutritive effects on healthy and varroa-parasitized bees. BMC Genomics. 2011; 12:496. https://doi.org/10.1186/1471-2164-12-496 PMID: 21985689

73. Amdam GV, Hartfelder K, Norberg K, Hagen A, Omholt SW. Altered physiology in worker honey bees (Hymenoptera: Apidae) infested with the mite Varroa destructor (Acari: Varroidae): a factor in colony loss during overwintering? J Econ Entomol. 2004; 97(3):741-7. PMID: 15279246

74. Asea A. Chaperokine-induced signal transduction pathways. Exerc Immunol Rev. 2003; 9:25-33. PMID: 14686091

75. Gorman MJ, Paskewitz SM. Serine proteases as mediators of mosquito immune responses. Insect biochemistry and molecular biology. 2001; 31(3):257-62. PMID: 11167095

76. Kanost MR, Jiang H, Wang Y, Yu XQ, Ma C, Zhu Y. Hemolymph proteinases in immune responses of Manduca sexta. Adv Exp Med Biol. 2001; 484:319-28. PMID: 11419000

77. Evans JD. Beepath: an ordered quantitative-PCR array for exploring honey bee immunity and disease. J Invertebr Pathol. 2006; 93(2):135-9. https://doi.org/10.1016/j.jip.2006.04.004 PMID: 16737710

78. Klaudiny J, Albert S, Bachanova K, Kopernicky J, Simuth J. Two structurally different defensin genes, one of them encoding a novel defensin isoform, are expressed in honeybee Apis mellifera. Insect Biochem. Mol. Biol. 2006; 35, 11-22.

79. Evans JD, Aronstein K, Chen YP, Hetru C, Imler J-L, Jiang M, Kanost M, Thompson GJ, Zou Z, Hultmark D. Immune pathways and defence mechanisms in honey bees Apis mellifera. Insect Mol Biol. 2006; 15:645-656. https://doi.org/10.1111/j.1365-2583.2006.00682.x PMID: 17069638

80. Amdam et al. Hormonal control of the yolk precursor vitellogenin regulates immune function and longevity in honeybees. Experimental Gerontology. 2004; 39: 767-772. https://doi.org/10.1016/j.exger.2004.02.010 PMID: 15130671

81. Yang X, Cox-Foster DL. Impact of an ectoparasite on the immunity and pathology of an invertebrate: evidence for host immunosuppression and viral amplification. Proc Natl Acad Sci U S A. 2005; 102 (21):7470-5. https://doi.org/10.1073/pnas.0501860102 PMID: 15897457

82. Steiner H. Peptidoglycan recognition proteins: on and off switches for innate immunity. Immunol Rev. 2004; 198:83-96. PMID: 15199956

83. Kaneko T, Silverman N. Bacterial recognition and signaling by the Drosophila IMD pathway. Cell Microbiol. 2005; 7:461-469. https://doi.org/10.1111/j.1462-5822.2005.00504.x PMID: 15760446

84. Sun W, Shen Y-H, Yang W-J, Cao Y-F, Xiang Z-H, Zhang Z. Expansion of the silkworm GMC oxidoreductase genes is associated with immunity. Insect biochemistry and molecular biology. 2012; 42 (12):935-45. https://doi.org/10.1016/j.ibmb.2012.09.006 PMID: 23022604

85. Brandt CR, Adang MJ, Spence KD. The peritrophic membrane: ultrastructural analysis and function as a mechanical barrier to microbial infection in Orgyia pseudotsugata. Journal of Invertebrate Pathology. 1978; 32(1):12-24.

86. Parker R, Guarna MM, Melathopoulos AP, Moon K-M, White R, Huxter E, et al. Correlation of proteome-wide changes with social immunity behaviors provides insight into resistance to the parasitic mite, Varroa destructor, in the honey bee (Apis mellifera). Genome Biol. 2012; 13(9):R81-2012. https://doi.org/10.1186/gb-2012-13-9-r81 PMID: 23021491

87. Foret S, Maleszka R. Function and evolution of a gene family encoding odorant binding-like proteins in a social insect, the honey bee (Apis mellifera). Genome Res. 2006; 16(11):1404-13. https://doi.org/10.1101/gr.5075706 PMID: 17065610

88. Zhang Y, Liu X, Zhang W, Han R. Differential gene expression of the honey bees Apis mellifera and A. cerana induced by Varroa destructor infection. J Insect Physiol. 2010; 56(9):1207-18. https://doi.org/10.1016/j.jinsphys.2010.03.019 PMID: 20346951