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. Amedian of 8.6 X 10⁵spores/ml was obtained.

Inoculation of bees with N. ceranaeand 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).

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 beebase 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.

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 days10 and 15 P.I. N. ceranae was not detected in any bees from the "control box". After day14 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 15694and 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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

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

1.3. Gene expression characteristics at days 5, 10 and 15 P.I. with N. ceranae.

Principle 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 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) and10 days P.I versus 15 days P.I (p = 1.92e-33).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

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

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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

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).

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 Mrjp6were 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.