Skip to main content


Differential gene expression and alternative splicing in insect immune specificity

Article metrics



Ecological studies routinely show genotype-genotype interactions between insects and their parasites. The mechanisms behind these interactions are not clearly understood. Using the bumblebee Bombus terrestris/trypanosome Crithidia bombi model system (two bumblebee colonies by two Crithidia strains), we have carried out a transcriptome-wide analysis of gene expression and alternative splicing in bees during C. bombi infection. We have performed four analyses, 1) comparing gene expression in infected and non-infected bees 24 hours after infection by Crithidia bombi, 2) comparing expression at 24 and 48 hours after C. bombi infection, 3) determining the differential gene expression associated with the bumblebee-Crithidia genotype-genotype interaction at 24 hours after infection and 4) determining the alternative splicing associated with the bumblebee-Crithidia genotype-genotype interaction at 24 hours post infection.


We found a large number of genes differentially regulated related to numerous canonical immune pathways. These genes include receptors, signaling pathways and effectors. We discovered a possible interaction between the peritrophic membrane and the insect immune system in defense against Crithidia. Most interestingly, we found differential expression and alternative splicing of immunoglobulin related genes (Dscam and Twitchin) are associated with the genotype-genotype interactions of the given bumblebee colony and Crithidia strain.


In this paper we have shown that the expression and alternative splicing of immune genes is associated with specific interactions between different host and parasite genotypes in this bumblebee/trypanosome model.


Invertebrate ecological studies have found infection outcomes within a given host-parasite system are variable. Parasite genotypes vary in their ability to infect, host genotypes vary in their ability to defend against infection, and infection outcome is based on the combination of both hosts’ and parasites’ genotypes [13]. That is, some interactions between invertebrate hosts and their parasites show specificity [4]. How is this level of specificity generated? An obvious answer would be an interaction between the parasite and the host’s immune response. This cannot be taken for granted however. Various ecological measures of disease outcome have been used to quantify genotype-genotype interactions. These measures include host mortality, fecundity and infection rate. Such measures cannot test directly if it is the immune response that produces this level of specificity [5]. It may be other non-immune processes, such as behavioural modification, reproductive isolation or life-history strategy changes, could explain such outcomes.

The interaction of the trypanosomal gut parasite Crithidia bombi with its host, the bumblebee, Bombus terrestris displays host x parasite genotype-genotype interactions [6, 7]. In this system, infection success depends on which strain is infecting which colony (representing very different genotypic backgrounds). There is evidence that the immune system has a role in generating this host-parasite specific response. A number of studies have found differential immune gene expression in response to Crithidia[810]. Previous work in our lab found increased Crithidia loads in bees whose expression of antimicrobial peptides was knocked down by RNAi [11]. It has even been shown that bees from different host genotypes induce differential expression of antimicrobial peptides (AMPs), according to the strain of C. bombi they had been infected with [12], that is specificity was found in the immune response itself. A recent paper using RNA-Seq found numerous genes are differentially expressed in a genotype-genotype fashion [13].

Although ideally we would have separate measures of infection levels, we could not do this due to a limited number of age controlled bees per colony. Based on their presumed importance in fighting Crithidia infections as mentioned above, we used antimicrobial peptide (AMP) expression as a proxy for a strong immune response in this manuscript.

Here, we carry out a transcriptome-wide analysis of gene expression in bees during C. bombi infection (two bumblebee colonies by two Crithidia strains). We have carried out four analyses, comparing 1) expression in infected and non-infected bees 24 hours after infection by Crithidia bombi (Infected versus uninfected) 2) expression at 24 and 48 hours after C. bombi infection (24 versus 48 hour), 3) determining differential gene expression associated with the host-parasite genotype-genotype interaction at 24 hours post infection (Specificity) and 4) determining alternative splicing associated with the host-parasite genotype-genotype interaction at 24 hours post infection (Specificity). Enrichment analysis was also carried out on differential expression data to see which categories of molecules are differentially regulated during infection. The results confirm our previous findings of up-regulation in antimicrobial peptide expression and provide a comprehensive overview of changes in and the specificity of gene expression and alternative splicing after exposure to 2 strains of C. bombi.

Results and discussion

The bumblebee colonies (host) and Crithidia bombi strains (parasite) used during this experiment are as previously described [12]. We have chosen samples that displayed a reciprocal pattern of expression for the three antimicrobial peptides (AMPs) tested in that paper. These were colony K5 (called K from now on) and Q1 (Q) and strains 6 and 8. K-8 showed a high expression in each of these AMPs, Q-8 a low expression level, Q-6 a high level and K-6 a low level of AMP expression. Bees were infected with various strains of Crithidia for either 24 or 48 hours. RNA samples were then collected and RNA-seq was carried out on these samples.

The sequences, statistics and annotations for all differentially expressed genes in each of the three differential expression analyses are available in Additional file 1: Table S1.

Genes differentially expressed at 24 hours post-infection (Infected versus uninfected)

Here we determined transcripts that are differentially expressed upon infection with Crithidia at 24 hours post-infection (Infected versus uninfected). We used both colonies and bees were infected with either strain 6 or 8. The statistical model was (0+colony+infect(yes/no)).31,843 unique transcripts were mapped to the transcriptome. 489 transcripts were found to be differentially expressed 24 hours post-infection (FDR < 0.05), including 324 downregulated and 165 upregulated transcripts. Reannotating the transcripts using Blast2GO (blastx against the nr database with e < 0.001), 109 had no BLAST hits. A further 68 had uninformative BLAST hits (anonymous predicted protein). The remaining 312 were used in the enrichment analysis. Figure 1 shows a summary of the enriched GO terms found (Fisher’s test p < 0.05). Defense response (GO:0006952, FDR = 0.047) and chitin metabolism (GO:0006030, FDR = 0.032) were the only processes significantly enriched at a more stringent level (FDR < 0.05).

Figure 1

Enriched GO terms found for differentially expressed genes at 24 hours post infection (infected versus uninfected). Using Blast2Go, we carried out an enrichment analysis (Fisher exact test p < 0.05) to see which GO terms are overrepresented relative to the entire transcriptome. These enriched GO terms were then summarized using Revigo.

Peritrophic membrane

The peritrophic matrix (PM) forms a layer composed of chitin and glycoproteins that lines the insect midgut lumen [14]. The PM facilitates digestion and forms a protective barrier to prevent the invasion of ingested pathogens [14, 15]. Fibrillin 1 (BTT14121_1), a venom protein precursor (BTT32193_1), Neurotrypsin (BTT07956_1), Peritrophin-1-like (BTT01709_1, BTT22959_1, BTT37215_1, BTT42262_1) and four chitinase transcripts (Chitinase 3: BTT23997_1 BTT38724_1, Chitinase 4 BTT20684_1, BTT23469_1) are downregulated upon infection. Fibrillins are extracellular matrix macromolecules, ubiquitous in the connective tissues [16]. BTT32193_1 was classed as a venom protein, but was also very similar to Chitinase 3 (blastx e = 1e-16). Chitinases modulate the structure and porosity of the PM [17]. Neurotrypsin is a serine protease expressed in the nervous system [18]. However in the protease domain it shares similarities with Sp22D, a chitin binding serine protease [19]. The chitin fibrils of the PM are assembled into a wide cross-hatched pattern connected by peritrophins [17]. A second group made up of Peritrophin-1 (BTT05886_1, BTT20661_1) and 3 further chitinase transcripts (Chitinase 2: BTT23246_1, Chitinase 3: BTT39163_1, Chitinase 4: BTT05313_1) is upregulated. Figure 2 shows the correlation of expression patterns between these sixteen transcripts related to chitin metabolism. There is some clustering, but not of any clear functional groups. Taken together this differential expression suggests an important role for the repair or restructuring of the peritrophic matrix in the bumblebees’ response to Crithidia.

Figure 2

Correlations of the chitin transcripts’ expression patterns that where differentially expressed twenty four hours post-infection compared to uninfected samples (infected versus uninfected). Clustering is produced based on Euclidean distances. The histogram shows the distribution of correlation values.


When the BLAST searches against the IIID and nr databases were combined, we found that 89 transcripts relate to canonical insect immune genes. We describe them in the order receptors, serine proteases, signalling pathways and effectors [4].

The Down syndrome cell adhesion molecule (Dscam), a pattern recognition receptor containing immunoglobulin domains has come to the forefront of research into insect immune specificity. Through alternative splicing it can generate thousands of different splice forms which would potentially allow an insect to recognise subtle differences in parasites. Its expression has been associated with infections in Drosophila and mosquitoes [20]. We found five downregulated transcripts annotated as immunoglobulin superfamily (Dscam included in hit list) (BTT03519_1, BTT08682_1, BTT15814_1, BTT26724_1, BTT27678_1) and one upregulated transcript (BTT03519_1).

Serine proteases

Serine proteases are important proteolytic enzymes in many molecular pathways. When these serine proteases are no longer needed, they are inactivated by serine protease inhibitors [21]. CLIP domain serine proteases mediate insect innate immunity [22]. Twenty one transcripts related to serine proteases, serine protease homologues or serine protease inhibitors were differentially expressed upon infection (see Table 1). Lipophorin receptor 2 (downregulated BTT34617_1) binds with serpins to aid in their endocytosis [23].

Table 1 List of transcripts associated with serine proteases and serine protease inhibitors found to be differentially expressed twenty four hours post infection with Crithidia bombi (Infected versus uninfected)

Signalling pathways

The insect immune system is regulated by three major signalling patways: Toll, Imd and JAK/STAT [24].

We found a transcript for Spatzle (BTT19738_1) downregulated at this time point. Activation of the Toll immune pathway requires the activation of Spatzle[24]. MyD88 (upregulated BTT15687_1) is a death domain-containing adaptor activated by Toll leading to the activation of Pelle. Dorsal (BTT25273_1) was also downregulated. The nuclear translocation of Dorsal, a member of the NF-kB family, in the Toll pathway induces the expression of many immune genes. We found an upregulated transcript (BTT09662_1) for Helicase89B part of the Toll and Imd Pathway. It is required downstream of NF-kB for the activation of AMP genes in Drosophila melanogaster[25]. ird5 (BTT03904_1 downregulated) codes for a catalytic subunit of an IkappaB kinase that cleaves Relish. Relish, a NF-KB factor in the Imd pathway, is an essential regulator of antimicrobial peptide gene induction.

In mammals semaphorins are crucially involved in various aspects of the immune response [26]. A semaphorin-5A-like transcript (BTT01850_1) was downregulated 24 hours post-infection. Semaphorin regulates the activity of Ras-family small GTPases [26]. A Ras-like protein11B transcript (BTT05368_1) was also downregulated. The Ras/MAPK pathway was found to be essential for the suppression of the Imd immune pathway in Drosophila [27].

Drumstick (downregulated BTT13062_1) interacts with the JAK/STAT pathway during its’ development role [28], but we could not find any information about its immune role. Two transcripts (BTT11590_1, BTT14205_1) of Puckered were downregulated. Puckered, which codes for a dual specificity phosphatase, is a key regulator of the c-Jun-N-terminal kinase (JNK) immune pathway [29]. Mpk2/p38a (downregulated BTT05769_1) is involved in the JNK Pathway and JAK/STAT Pathway. Heat-shock factor activation by p38 is a recently discovered part of antimicrobial reactions in Drosophila[30]. We found two heat shock protein transcripts (BTT23758_2, BTT37030_1) and one other (BTT17701_1) that were downregulated and upregulated respectively. These are all involved in the JAK/STAT pathway.


In our previous paper [8] we found that antimicrobial peptides were upregulated at 24 hours post-infection. We would expect the same trend here. Indeed, we found that five transcripts for defensin (BTT06274_2, BTT8490_1, BTT10405_1, BTT14019_1, and BTT42034_1) and three hymenoptaecin transcripts (BTT18071_1, BTT24170_1, BTT24170_2) were all upregulated. An apidaecin precursor (BTT33652_1) was downregulated. Apidaecin has recently been shown to be expressed in bumblebees [31] including in response to Crithidia[9, 10, 13, 32]. The downregulated beta-amyloid-like protein (BTT20240_1) has been shown to be an antimicrobial peptide in mammals [33]. Hemolectin (BTT15326_1, upregulated) is a clotting protein known to have a role against gram-negative bacteria [34].

Reactive oxygen species (ROS) are generated by respiration in the mitochondria or as part of the immune response [35]. P450 cytochromes are oxidases acting terminally in monooxygenase systems [36]. Some are regulated in response to infection possibly either as direct immune responders [37], producing nitric oxide (NO) or other reactive oxygen radicals or as part of the host detoxification process decreasing oxidative stress after an infection [35]. A number of cytochromes P450 were differentially expressed 24 hours post infection. Ten cytochrome p450 transcripts (Cyp4p3: BTT05294_1, BTT20848_1, BTT22253_1, BTT23317_1, BTT32674_1, cytochrome P450 4g15: BTT23811_1, BTT32459_1, cytochrome P450 6k1: BTT35547_1, BTT40653_1, cytochrome P450 6a14: BTT38445_1) were found to be downregulated. Three other cytochrome P450 transcripts (Cyp4p3: BTT21216_1, BTT35543_1, cytochrome P450 315a1: BTT26726_1) were upregulated. Several other cytochromes (cytochrome b: BTT20524_1, BTT39776_1, BTT41896_1, and cytochrome c: BTT05255_2) were downregulated.

Numerous other actors in the production of ROS were found to be differentially expressed. TPX4 (BTT13285_1), coding for a Thioredoxin-dependent peroxidase, which detoxifies hydrogen peroxide was downregulated. This gene was found be differentially expressed during Plasmodium infection in Anopheles gambiae[38]. Calcineurin (BTT08150_1, BTT26273_1) was found to be downregulated 24 hours post-infection, which agrees with our previous findings [8]. In infected D. melanogaster larvae, NO signals are enhanced by Calcineurin to promote induction of systemic immune responses via the Imd signalling pathway [39].

We found downregulation of sortilin-related receptor-like (BTT31654_1). In mammals, sortilin aids in phagocytosis [40]. Two downregulated transcripts (BTT35021_1, BTT08756_1) were matched to croquemort, which codes for a key scavenger receptor although of apoptotic cells rather than parasites [41]. Annexin IX (downregulated BTT02025_1) has been shown to be induced by septic injury in Drosophila and is thought to encode for an anticoagulant [42].


Major royal jelly protein (BTT05317_2, BTT36365_1 upregulated) has been shown to have antimicrobial properties and to be expressed in response to bacterial infection in honeybees [43, 44]. Vitellogenin (downregulated BTT36006_1) is a potent regulator of the immune response in honeybees [45]. Several orthologs of putative Drosophila immune loci were differentially expressed 24 hours post-infection (CG12505: BTT00934_1, CG18348: BTT04397_1, CG7296: BTT15035_1, BTT18395_1, CG8791: BTT18908_1, CG5527: BTT35653_1, Fst: BTT11511_1). The downregulated CG4393 (BTT05817_1) is weakly analogous to TNF receptor associated factor 3 (TRAF3) that mediates signal transduction involved in mammalian immune responses. Downregulated BTT37289_1 codes for a putative fatty acyl-CoA reductase.

Genes differentially expressed between 24 hours post-infection and 48 hours post-infection (24 versus 48 hours)

Here we determined the gene expression difference between 24 hours and 48 hours post infection (24 versus 48 hours). Both colonies were used but only strain six was used for infection. The statistical model was (0+colony + time).

43 transcripts were differentially expressed between 24 hours post-infection and 48 hours post-infection. Of these 17 had no BLAST hits. A further six had uninformative BLAST hits (anonymous predicted protein). The remaining 20 were used in the analysis. Defence response was the only GO term significantly enriched compared to the whole transcriptome (Fisher exact test: FDR = 0.00015), with seven transcripts. Three transcripts correspond to Hymenoptaecin (BTT18071_1, BTT24170_1, BTT24170_2). They were all upregulated. This suggests a continuing strong AMP production 48 hours after infection. This agrees with other immune assays in bumblebees [46]. Argonaute-2, a RNA-silencing endonuclease, is involved in antiviral defence in insects (downregulated BTT02484_1) [47]. GstD8, a glutathione S-transferase, is involved in the detoxification process (upregulated BTT04810_1) [48]. Dopa decarboxylase (upregulated BTT28048_1) converts L-dopa to dopamine during the melanisation process [49]. SCR-B9 (upregulated BTT40924_1) codes for a scavenger receptor protein. Scavenger receptor proteins have been found to be microbial pattern recognition receptors in Drosophila[50].

Genes differentially expressed depending on host genotype – parasite genotype interactions (Specificity)

We choose samples that displayed a reciprocal pattern of expression for the three antimicrobial peptides (AMPs) (Defensin, Abaecin and Hymenoptaecin) in a previous study [12]. These were colony K and Q and strains 6 and 8. K-8 showed a high AMP expression, Q-8 a low expression level, Q-6 a high level and K-6 a low level of AMP expression. All data was collected at 24 hours post infection. The statistical model was (0+colony*strain).

There were 591 differentially expressed transcripts (FDR < 0.05). Reannotating the transcripts using Blast2GO (blastx against the nr database with e < 0.001), 150 had no BLAST hits. A further 64 had uninformative BLAST hits (anonymous predicted protein). There were 109 transcripts that had previously been found to be differentially expressed at 24 hours post infection.

Of the 591 transcripts, 132 were upregulated and 459 were downregulated. Up or downregulation does not have the same meaning here as in the infected versus uninfected model were there was a clear baseline (uninfected). Our model used colony K strain 8 as the final contrast. From our previously published qPCR data [12], we know the colony K strain 8 interaction displayed the highest levels of AMPs (effectors). Therefore when we say a transcript is upregulated, we mean it is upregulated in this high immune response interaction.

As with the infection data, we combined the BLAST searches against the IIID and nr databases. Ninety transcripts correspond to canonical insect immune genes. We again describe them in the order receptors, serine proteases, signalling pathways and effectors.


Two transcripts were associated with Gram-negative binding proteins (upregulated GNBP, BTT03533_1 and downregulated GNBP1-2 BTT35513_1). Although GNBPs are most associated with defense against gram-negative bacteria, they have been show to have a role in response to Plasmodium infections [51]. C-type lectins (CTLs) bind carbohydrates and mediate processes including pathogen recognition [52]. CTL4 is agonist to Plasmodium infections in mosquitoes [52]. A CTL4 transcript (BTT29328_1) was found to be downregulated.

One downregulated transcript was related to Dscam (BTT12755_1). A further fourteen downregulated transcripts were part of the Ig superfamily (IGFn3-1: BTT05561_1, BTT05581_1, BTT08682_1, BTT12655_1, BTT13442_1, BTT14516_1, BTT18750_1, BTT21156_1, BTT22598_1, BTT22819_1, BTT23339_1, BTT24070_1, IGFn3-7: BTT08109_1, BTT09498_1) and one was upregulated (IGFn3-8: BTT03519_1). Dscam and most of the other Ig superfamily transcripts cluster together in the top right of Figure 3, suggesting they are similarly expressed.

Figure 3

Correlations of the immunity transcripts’ expression patterns that were differentially expressed depending on host genotype – parasite genotype interactions (Specifcity). Clustering is produced based on Euclidean distances. The histogram shows the distribution of correlation values.

Serine proteases

28 transcripts related to serine proteases, serine protease homologues or serine protease inhibitors were differentially expressed (see Table 2).

Table 2 List of transcripts associated with serine proteases and serine protease inhibitors found to be differentially expressed in the specificity analysis

Signalling pathways

The Toll-like receptor 18Wheeler (BTT35732_1) and Toll 10 (BTT09386_1) were both upregulated. 18Wheeler has been shown to be important in the anti gram-negative immune response in Drosophila larvae [53]. Dorsal 1A (BTT04010_1), a transcription factor that is an important part of the Toll pathway, was downregulated. A transcript for Spatzle 1-2 was also downregulated (BTT10679_1).

The tyrosine kinase Pvr (BTT04822_1), which inhibits JNK activation [54] was downregulated. Jun, a transcription factor of the JNK pathway was downregulated (BTT13636_1). Mpk2/p38a (downregulated BTT16580_1) and MAPKKK9 (downregulated BTT04404_1) are mitogen-activated protein kinases involved in the JNK Pathway and JAK/STAT pathways. We found two heat shock protein transcripts (BTT17371_1, BTT22195_1) and one other (BTT17701_1) that were downregulated and upregulated respectively. These are all involved in the JAK/STAT pathway. Akt 1 (downregulated BTT14188_1) is part of the insulin/insulin-like growth factor 1 signaling (IIS) cascade. IIS plays a critical role in the regulation of innate immunity. Activation of Akt signaling leads to a decrease in malaria infection intensity in mosquitoes [55].


Five transcripts relate to the AMPs defensin (BTT06274_2, BTT42034_1) and hymenoptaecin (BTT18071_1, BTT24170_1, BTT24170_2). They were all upregulated. An apidaecin precursor (BTT20828_1) was upregulated. Hemolectin had three downregulated transcripts (BTT14194_1, BTT17013_1, BTT26614_1) and one upregulated (BTT15326_1). Argonaute-2, a RNA-silencing endonuclease, is involved in antiviral defense in insects (downregulated BTT02374_1) [47].

Eater encodes a transmembrane receptor involved in phagocytosis in Drosophila[56]. A transcript (BTT11132_1) relating to Eater was upregulated. The melanisation process component Dopa decarboxylase (BTT19093_1) was upregulated. Another enzyme involved in melanisation, laccase was found to be downregulated (BTT20241_1, BTT33633_1) [57].

A cytochrome P450 transcript Cyp4p3 BTT40653_1 was upregulated. Two previously unseen Cyp4p3 transcripts (BTT05254_1, BTT20622_2) were upregulated and one (BTT36257_1) downregulated. TPX4 (BTT13285_1) that codes for a Thioredoxin-dependent peroxidase was downregulated.


A small number of transcripts were related to chitin metabolism. SCRASP1 has a chitin-binding domain that has been hypothesized to sense chitin in response to injury and to transduce signals via the serine protease domain [58]. We found an upregulated transcript related to SCRASP 1 (BTT41923_1). A peritrophin precursor was also upregulated (BTT10727_1), as was a chitinase 3 transcript (BTT23246_1).

Retinoid and fatty-acid-binding protein (RfaBp) (BTT07678_1) was downregulated. RfaBp was found to be upregulated upon injection of LPS in Drosophila during a proteomic study [59]. Notch (upregulated BTT09545_1) is involved is involved in many signaling events in development including the specification of crystal cells in Drosophila melanogaster[60]. Finally, several orthologs of putative Drosophila immune loci were found to be differentially expressed (CG5527: BTT08512_1, CG12505: BTT00934_1, CG13323: BTT38025_1, BTT38087_1, CG17560: BTT02877_1 downregulated, BTT05845_1CG18348: BTT20843_1).

Genes alternatively spliced depending on host genotype – parasite genotype interactions (Specificity)

Using an interaction model in DEXSeq we identified genes which show differential exon usage due to the interaction of strain and colony (FDR < 0.01). The complete output of the DEXSeq analysis is available in the Additional file 2: Table S2. 615 loci displayed alternative splicing depending on the the interaction between the host genotype and the parasite genotype. The sequences, statistics and annotations for these loci is available in the Additional file 3: Table S3. 98 of the loci had a significant match against the Insect Innate Immunity Database (IIID). Eleven of these are related to receptors including lectins (XLOC_007616 XLOC_007614 XLOC_007613, XLOC_012985, XLOC_011830 XLOC_011825), beta-1,3-glucan recognition protein (XLOC_003146 XLOC_003142 XLOC_003143) and seven transcripts relating to immunoglobulin and fibronectin domains (XLOC_006019, XLOC_006020, XLOC_010845 XLOC_010849, XLOC_012355, XLOC_013957, XLOC_012379, XLOC_014287). Several of these transcripts containing immunogloblin and fibronectin domains returned Dscam as a blast hit. The transcript with the largest number of alternatively spliced exons was XLOC_010845 XLOC_010849 with 25 variable exons out of 71 total exons (see Figure 4). Against the Bombus terrestris genome a blast search of XLOC_010845 XLOC_010849 returned a twitchin-like gene (XM_003393691). XLOC_010845 XL OC_010849 aligns with the first 31,160 base pairs of the twitchin-like gene (80,078 bp in total). Our cufflinks Twitchin gene model matches almost perfectly the one generated automatically by the Bombus genome team.

Figure 4

Twitchin (XLOC_010845 XLOC_010849) differential exon usage for each of the four host-parasite strain combinations (K6, K8 Q6, Q8). DEXSeq removes the gene-level changes in expression so the differential exon usage is clear. Below is the gene model produced by our cufflinks analysis. The grey blocks are normal exons, the purple blocks represent those exons showing alternative splicing.


We present a comprehensive transcriptomic analysis of gene expression in this important model host-parasite system. We have identified 489 bumblebee genes whose expression are changed upon infection with Crithidia. We also found 591 genes whose expression is associated with the interaction between host and parasite genotypes and therefore show specificity in their expression patterns. Six hundred and fifteen genes showed alternative splicing in response to parasite-strain interactions.

Our AMP expression data is consistent with the proposed importance of antimicrobial peptides in the specific defence against Crithidia[8, 11, 12]. It is also clear that several other effectors including ROS and phagocytosis may be important. Several immune pathways seem to be invovled in the anti-Crithidia response. These include the Toll, Imd and JAK/STAT pathways.

There are a larger proportion of receptor transcripts showing differential expression found in the specificity analysis (3.2% 19/591) compared to the infection analysis (1.2% 6/489). This is not surprising, as it is be expected that a specific immune response to a given strain would be based mainly on how it is recognised. Although several receptors, including GNBPs and lectins, are differentially expressed, the most exciting discovery is the large number of transcripts related to Dscam. The Down syndrome cell adhesion molecule (Dscam), a pattern recognition receptor has come to the forefront of research into insect immune specificity as thousands of different splice forms are generated and it is associated with insect immunity [20]. In the fruit fly Drosophila, silencing of Dscam retards the insect’s capacity to engulf bacteria by phagocytosis [61]. In mosquito Anopheles, the Dscam splice forms produced in response to parasite exposure differs between bacteria and the malarial causitative agent Plasmodium and between Plasmodium berghei and Plasmodium falciparum[62]. This has been tempered by the finding that Dscam diversity does not increase with exposure to increasing heterogeneity of Plasmodium falciparum genotypes [20]. Recently it has been shown that Dscam specificity is mediated by the transcriptional regulation of specific splicing factors downstream of the activation of the Toll and IMD pathways [63]. Our results suggest that Dscam related genes may be important in differentiating strains of the trypanosome Crithidia bombi.

The alternative splice analysis also found a number of receptor genes. This included numerous Dscam related genes. This is encouraging as alternative splicing is the mechanism through which Dscam generates the variation that is thought to be useful for immune recognition [20]. The gene with the largest number of alternatively spliced exons was Twitchin. This gene was also downregulated 24 hours post- infection (BTT27678_1). Five different transcripts of Twitchin (BTT12655_1, BTT13442_1, BTT21156_1, BTT22598_1, BTT23339_1) were expressed in a genotype-genotype fashion in the specificity differential expression analysis. Like Dscam, Twitchin possesses a large nymber of fibronectin and immunoglobulin domains. Twitchin is part of the titin family of genes. They produce large filamentous proteins that mediate the transduction of mechanical signals in muscles [64]. However a titin gene was found to show differential exon usage depending on if a mosquito was infected with bacteria or Plasmodium[65]. Twitchin is an exciting possible candidate gene for the source of specifcity in this system and could be a fruitful avenue of research.

A large number of genes are downregulated after infection. Naively, we might expect genes fighting infections to be upregulated after infection. We have two possible explanations. The more exciting of the two is that the parasites are modifying the immune response. Trypanosoma cruzi, the causative agent of Chagas disease has been shown to do this in it’s insect host [66]. The other is that for a lot of the earlier responder genes such as receptors we may have missed the peak of their expression with our earliest timepoint being at twenty four hours post-infection. Repeating at earlier timepoints would establish if that was the case.

We infected the bees with faceces from other bees. We chose this method over in vitro culturing to prevent possible attentuation of strains’ infectivity associated with culturing [67]. One possibility is that by using faeces we may be introducing hidden infections or gut microbiota from the donor queens. We have attempted to mitigate this by using bees from a single colony to culture and grow the Crithidia. Although the queens faeces may indeed contain hidden infections and microbiota, they all must be passed through the same host background before they are used experimentally.

We found a number of genes associated with chitin metabolism differentially regulated 24 hours post-infection. Through several pathways chitin metabolism is fundamental to invertebrate immunity [68]. As an aside, an intriguing hypothesis is that chitin metabolism is the nexus through which defense against predators and against parasites are traded-off [68]. Our data suggests that the peritrophic matrix may be fundamental in the bee’s defence against Crithidia. The peritrophic matrix acts as an immunological barrier against trypanosomes. Tsetse flies with an underdeveloped PM have lower levels of refractoriness to trypanosome infections [69]. This is due to a premature immune response; the trypanosomes get through the PM quicker and stimulate the immune response at an earlier stage compared to refractory Tsetse flies.

A recently published paper by one of the authors, found genotype x genotype interactions in the expression of a smaller number of genes [13]. We hypothesise that our much larger catalogue of genes, including Dscam and Twitchin, is due to our different experimental design. They had four colonies and three parasite strains all from queens caught in the same local area. We choose samples that displayed a 2 (colony) × 2 (parasite strain) reciprocal pattern of expression for the three antimicrobial peptides (AMPs) (Defensin, Abaecin and Hymenoptaecin) in a previous study [12]. These were colony K and Q and strains 6 and 8. K-8 showed a high AMP expression, Q-8 a low expression level, Q-6 a high level and K-6 a low level of AMP expression. The two colonies used were from different populations, one wild caught and one commercial. This increased the potential differences in their response to the two strains. In turn, this increased our likelihood of detecting differential expression and exon usage using RNA-seq.

In this paper we have shown that the expression and alternative splicing of immune genes is associated with specific interactions between different host and parasite genotypes in this bumblebee/trypanosome model. Future RNAi work could knockdown candidate genes thereby altering these specific interactions to directly examine their biological significance.


Sample collection

The bumblebee colonies (host) and Crithidia bombi strains (parasite) used during this experiment are as previously described [12]. We have chosen samples that displayed a reciprocal pattern of expression for the three antimicrobial peptides (AMPs) tested in that paper. These were colony K5 (called K from now on) and Q1 (Q) and strains 6 and 8. K-8 showed a high expression in each of these AMPs, Q-8 a low expression level, Q-6 a high level and K-6 a low level of AMP expression. Experiments were carried out on one commercially reared bumblebee colony from Koppert Biological Systems U.K. (Colony K) and one colony from a wild caught queen (Colony Q). Faecal samples from these colonies were checked under a light microscope to ensure there was no Crithidia bombi present [70]. All parasite isolates used originated from wild queens collected in Spring 2008 in the University of Leicester botanic garden. The Crithidia from each individual queen was infected into a group of 10 workers from a different colony to amplify the strain and to provide a source for experimental infections. Experiments began when the colonies had a minimum of thirty workers, approximately four weeks old. Between observations, colonies were fed ad libitum with pollen (Percie du sert, France) and 50% diluted glucose/fructose mix (Meliose – Roquette, France). Before and during the experiments colonies were kept at 26°C and 60% humidity in constant red light.


To prepare C. bombi isolates, faeces was collected from infected workers and mixed with 50% diluted Meliose to create a standardized dose of 500 Crithidia cells per microlitre of inoculum. Previous studies had shown that such inocula, prepared from different colonies, are genotypically different [7] and generate specific responses in novel hosts [6]. We infected a sample of workers from each of K and Q bumblebee colonies (representing different host lines) with an inoculum of faeces from each of the two wild infected queens (6 and 8 Crithidia strain). We also collected uninfected controls, which were sacrificied at five days old (fours days plus 24 hours). Bees were four days old at the time of infection. Bees were collected over several days and distributed across treatment groups [71]. After infection bees were kept in colony × strain groups (1–3 individuals depending on day collected) and fed ad libitum. Twenty four hours or 48 hours post infection the bees were sacrificed by freezing in liquid nitrogen and stored at minus 80°C.

RNA sample preparation and sequencing

Total RNA was extracted from 23 individual homogenised abdomens using Tri-reagent (Sigma-Aldrich, UK). Samples (Colony-Strain-Timepoint (number of replicates)) were K-6-24 (3), K-6-48 (3), K-8-24 (3), K-8-48 (3), K-Uninfected (2), Q-6-24 (3), Q-6-48 (3), Q-8-24 (2), Q-uninfected (1). Any residual contaminants were removed from the RNA using the RNeasy mini kit (Qiagen, UK) and manufacturer’s RNA clean-up protocol. To remove residual genomic DNA, RNA samples were treated with DNase (Sigma-Aldrich, UK). TruSeq RNA-seq libraries were made from the 23 samples at NBAF Edinburgh. Sequencing was performed on an Illumina HiSeqxxxxxxxxR2000 instrument (Illumina, Inc.) by the manufacturer’s protocol. Multiplexed 50 base single-read runs were carried out yielding an average of 12M reads per sample.

Differential gene expression analysis

The reference transcriptome was downloaded from terrestris/ [31]. Functional annotation related to the transcriptome was obtained using the BLAST2GO package [72]. Alignment was done using GSNAP (version 2012-07-20) [73]. Only reads that mapped uniquely were selected for further analysis. Counts were generated per transcript for each sample.

Differential expression analysis was performed using the edgeR (3.4.0) package [74] in R (3.0.1) [75]. Normalization factors were computed using the TMM technique, after which tagwise dispersions were calculated and subjected to a generalized linear model (GLM). Resulting p values were subjected to Benjamini–Hochberg multiple testing correction to derive FDRs; only transcripts with a FDR < 0.05 were considered for further analysis. Three separate GLMs were carried out. One looked for transcripts that are differentially expressed upon infection with Crithidia at 24 hours post-infection (Infected versus uninfected) (0+colony+infect(yes/no)). “Infect” here are bees infected with either strain 6 or 8. Another GLM looked at the gene expression difference between 24 hours and 48 hours post strain 6 infection (24 versus 48 hours)(0+colony + time). The third GLM looked for transcripts that were expressed in a specific pattern at 24 hours post-infection (specifcity)(0+colony*strain).

Using Blast2Go, we then carried out an enrichment analysis (Fisher exact test) on each of these lists of differentially expressed genes to see which GO terms are overrepresented relative to the entire transcriptome. We then used REVIGO to summarize and visualise these terms [76].

For each of the lists of differentially expressed transcripts we also carried out a blastx analysis against the insect innate immunity database (IIID) [77]. We used the BLOSUM62 matrix with a word size of 3. The results were filtered to only contain hits with an E-value < 1e-10 and a bit score ≥30.

Alternative splice analysis

The eleven samples used in the specificity analysis above were also tested for alternative splicing.

Alignment and creation of gene set

Reads were first aligned to the Bombus terrestris reference genome (AELG00000000.1) using the fast splice junction mapper Tophat [78]. Preliminary sequence data was obtained from Baylor College of Medicine Human Genome Sequencing Center website at The Tophat produced alignment files were then passed to Cufflinks to generate a transcriptome assembly for each sample [79]. These assemblies were then merged together using the Cuffmerge utility [79], into a general transcriptome assembly.

DEXSEQ analysis

The BAM files from the Tophat analysis were converted into SAM format using SAMTools [80]. The GTF file from Cuffmerge was flattened into a GFF file with collapsed exon counting bins using the Python script found in HTSeq [81]. For each SAM file, (HTSeq) counts the number of reads that overlap with each of the exon counting bins defined in the flattened GFF file. We tested for differential exon usage using the R package DEXSeq [82]. We used a full model of counts = sample+exon+(colony+strain)*exon+(colony:strain)*I(exon==exonID) and a reduced model of counts = sample+exon+(colony+strain)*exon. This identified genes which show differential exon usage due to the interaction of strain and colony (FDR < 0.01).

Blast analysis

We extracted the nucleotide sequence for all differentially expressed transcripts and searched for any matching sequence on NCBI using BLASTn [83] with an E-value cutoff of 0.001, restricting the sequences to those from B. terrestris.

The protocol reported here conforms to the regulatory requirements for animal experimentation in the United Kingdom.


  1. 1.

    Carius HJ, Little TJ, Ebert D: Genetic variation in a host-parasite association: Potential for coevolution and frequency-dependent selection. Evolution. 2001, 5: 1136-1145.

  2. 2.

    Schmid-Hempel P: On the evolutionary ecology of host-parasite interactions: addressing the question with regard to bumblebees and their parasites. Naturwissenschaften. 2001, 88: 147-158. 10.1007/s001140100222.

  3. 3.

    Lambrechts L, Scott TW, Gubler DJ: Consequences of the expanding global distribution of aedes albopictus for dengue virus transmission. PLoS Negl Trop Dis. 2010, 4: 5-

  4. 4.

    Schmid-Hempel P: Natural insect host-parasite systems show immune priming and specificity: puzzles to be solved. Bioessays. 2005, 27: 1026-1034. 10.1002/bies.20282.

  5. 5.

    Hauton C, Smith VJ: Adaptive immunity in invertebrates: a straw house without a mechanistic foundation. Bioessays. 2007, 29: 1138-1146. 10.1002/bies.20650.

  6. 6.

    Schmid-Hempel P, Puhr K, Kruger N, Reber C, Schmid-Hempel R: Dynamic and genetic consequences of variation in horizontal transmission for a microparasitic infection. Evolution. 1999, 3: 426-434.

  7. 7.

    Schmid-Hempel P, Reber-Funk C: The distribution of genotypes of the trypanosome parasite, crithidia bombi, in populations of its host bombus terrestris. Parasitology. 2004, 129: 147-158. 10.1017/S0031182004005542.

  8. 8.

    Riddell CE, Sumner S, Adams S, Mallon EB: Pathways to immunity: temporal dynamics of the bumblebee (bombus terrestris) immune response against a trypanosomal gut parasite. Insect Mol Biol. 2011, 20 (4): 529-540. 10.1111/j.1365-2583.2011.01084.x. doi:10.1111/j.1365-2583.2011.01084.x,

  9. 9.

    Brunner FS, Schmid-Hempel P, Barribeau SM: Immune gene expression in bombus terrestris: Signatures of infection despite strong variation among populations, colonies, and sister workers. PLoS ONE. 2013, 8 (7): doi:10.1371/journal.pone.0068181. Accessed 2013-10-15,

  10. 10.

    Barribeau SM, Schmid-Hempel P: Qualitatively different immune response of the bumblebee host, bombus terrestris, to infection by different genotypes of the trypanosome gut parasite, crithidia bombi. Infect Genet Evol. 2013, 20: 249-256. doi:10.1016/j.meegid.2013.09.014. Accessed 2014-03-11,

  11. 11.

    Deshwal S, Mallon EB: Antimicrobial peptides play a functional role in bumblebee anti-trypanosome defense. Dev Comparative Immunol. 2014, 42 (2): 240-243. 10.1016/j.dci.2013.09.004. doi:10.1016/j.dci.2013.09.004. Accessed 2013-10-22,

  12. 12.

    Riddell C, Adams S, Schmid-Hempel P, Mallon EB: Differential expression of immune defences is associated with specific host-parasite interactions in insects. PLoS ONE. 2009, 4: 7621-10.1371/journal.pone.0007621.

  13. 13.

    Barribeau SM, Sadd BM, Schmid-Hempel P, Plessis L d: Gene expression differences underlying genotype-by-genotype specificity in a host–parasite system. Proc Natl Acad Sci U S A. 2014, 111 (9): 3496-3501. 10.1073/pnas.1318628111. doi:10.1073/pnas.1318628111. Accessed 2014-03-06,

  14. 14.

    Kuraishi T, Binggeli O, Opota O, Buchon N, Lemaitre B: Genetic evidence for a protective role of the peritrophic matrix against intestinal bacterial infection in drosophila melanogaster. Proc Natl Acad Sci U S A. 2011, 108 (38): 15966-15971. 10.1073/pnas.1105994108. doi:10.1073/pnas.1105994108,

  15. 15.

    Lehane MJ: Peritrophic matrix structure and function. Annu Rev Entomol. 1997, 42: 525-550. 10.1146/annurev.ento.42.1.525. doi:10.1146/annurev.ento.42.1.525,

  16. 16.

    Isogai Z, Ono RN, Ushiro S, Keene DR, Chen Y, Mazzieri R, Charbonneau NL, Reinhardt DP, Rifkin DB, Sakai LY: Latent transforming growth factor beta-binding protein 1 interacts with fibrillin and is a microfibril-associated protein. J Biol Chem. 2003, 278 (4): 2750-2757. 10.1074/jbc.M209256200. doi:10.1074/jbc.M209256200,

  17. 17.

    Dinglasan RR, Devenport M, Florens L, Johnson JR, McHugh CA, Donnelly-Doman M, Carucci DJ, Yates IIIJR, Jacobs-Lorena M: The anopheles gambiae adult midgut peritrophic matrix proteome. Insect Biochem Mol Biol. 2009, 39 (2): 125-134. 10.1016/j.ibmb.2008.10.010. doi:10.1016/j.ibmb.2008.10.010,

  18. 18.

    Gschwend TP, Krueger SR, Kozlov SV, Wolfer DP, Sonderegger P: Neurotrypsin, a novel multidomain serine protease expressed in the nervous system. Mol Cellular Neurosci. 1997, 9 (3): 207-219. 10.1006/mcne.1997.0616. doi:10.1006/mcne.1997.0616. Accessed 2013-10-17,

  19. 19.

    Danielli A, Loukeris TG, Lagueux M, Richman A, Müller H-M: A modular chitin-binding protease associated with hemocytes and hemolymph in the mosquito anopheles gambiae. Proc Natl Acad Sci U S A. 2000, 97 (13): 7136-7141. 10.1073/pnas.97.13.7136. doi:10.1073/pnas.97.13.7136. Accessed 2013-10-17,

  20. 20.

    Smith PH, Mwangi JM, Afrane YA, Yan G, Obbard DJ, Ranford-Cartwright LC, Little TJ: Alternative splicing of the anopheles gambiae dscam gene in diverse plasmodium falciparum infections. Malar J. 2011, 10 (1): 156-10.1186/1475-2875-10-156. doi:10.1186/1475-2875-10-156. Accessed 2013-04-19,

  21. 21.

    Zhao P, Dong Z, Duan J, Wang G, Wang L, Li Y, Xiang Z, Xia Q: Genome-wide identification and immune response analysis of serine protease inhibitor genes in the silkworm, bombyx mori. PLoS ONE. 2012, 7 (2): 31168-10.1371/journal.pone.0031168. doi:10.1371/journal.pone.0031168. Accessed 2013-11-06,

  22. 22.

    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 (5): 603-614. 10.1111/j.1365-2583.2006.00684.x. doi:10.1111/j.1365-2583.2006.00684.x. Accessed 2013-04-19,

  23. 23.

    Soukup SF, Culi J, Gubb D: Uptake of the necrotic serpin in drosophila melanogaster via the lipophorin receptor-1. PLoS Genet. 2009, 5 (6): 1000532-10.1371/journal.pgen.1000532. doi:10.1371/journal.pgen.1000532. Accessed 2013-10-08,

  24. 24.

    Lemaitre B, Hoffmann J: The host defense of drosophila melanogaster. Ann Rev Immunol. 2007, 25: 697-743. 10.1146/annurev.immunol.25.022106.141615.

  25. 25.

    Yagi Y, Ip YT: Helicase89b is a mot1p/BTAF1 homologue that mediates an antimicrobial response in drosophila. EMBO Rep. 2005, 6 (11): 1088-1094. 10.1038/sj.embor.7400542. doi:10.1038/sj.embor.7400542,

  26. 26.

    Takamatsu H, Kumanogoh A: Diverse roles for semaphorin plexin signaling in the immune system. Trends Immunol. 2012, 33 (3): 127-135. 10.1016/ doi:10.1016/ Accessed 2013-04-16,

  27. 27.

    Ragab A, Buechling T, Gesellchen V, Spirohn K, Boettcher A-L, Boutros M: Drosophila ras/MAPK signalling regulates innate immune responses in immune and intestinal stem cells. EMBO J. 2011, 30 (6): 1123-1136. 10.1038/emboj.2011.4. doi:10.1038/emboj.2011.4. Accessed 2013-04-16,

  28. 28.

    Johansen KA, Iwaki DD, Lengyel JA: Localized JAK/STAT signaling is required for oriented cell rearrangement in a tubular epithelium. Development. 2003, 130 (1): 135-145. 10.1242/dev.00202. doi:10.1242/dev.00202. Accessed 2013-04-16,

  29. 29.

    Karkali K, Panayotou G: The drosophila DUSP puckered is phosphorylated by JNK and p38 in response to arsenite-induced oxidative stress. Biochem Biophys Res Commun. 2012, 418 (2): 301-306. 10.1016/j.bbrc.2012.01.015. doi:10.1016/j.bbrc.2012.01.015. Accessed 2013-04-19,

  30. 30.

    Chen J, Xie C, Tian L, Hong L, Wu X, Han J: Participation of the p38 pathway in drosophila host defense against pathogenic bacteria and fungi. Proc Natl Acad Sci U S A. 2010, 107 (48): 20774-20779. 10.1073/pnas.1009223107. doi:10.1073/pnas.1009223107. Accessed 2013-04-19,

  31. 31.

    Colgan TJ, Carolan JC, Bridgett SJ, Sumner S, Blaxter ML, Brown MJ: Polyphenism in social insects: insights from a transcriptome-wide analysis of gene expression in the life stages of the key pollinator, bombus terrestris. BMC Genomics. 2011, 12 (1): 623-10.1186/1471-2164-12-623. doi:10.1186/1471-2164-12-623. Accessed 2013-04-05,

  32. 32.

    Brunner FS, Schmid-Hempel P, Barribeau SM: Protein-poor diet reduces host-specific immune gene expression in bombus terrestris. Proc Biol Sci /R Soc. 2014, 281 (1786): doi:10.1098/rspb.2014.0128,

  33. 33.

    Soscia SJ, Kirby JE, Washicosky KJ, Tucker SM, Ingelsson M, Hyman B, Burton MA, Goldstein LE, Duong S, Tanzi RE, Moir RD: The alzheimer’s disease-associated amyloid beta protein is an antimicrobial peptide. PLoS ONE. 2010, 5 (3): 9505-10.1371/journal.pone.0009505. doi:10.1371/journal.pone.0009505. Accessed 2013-04-16,

  34. 34.

    Lesch C, Goto A, Lindgren M, Bidla G, Dushay MS, Theopold U: A role for hemolectin in coagulation and immunity in drosophila melanogaster. Dev Comp Immunol. 2007, 31 (12): 1255-1263. 10.1016/j.dci.2007.03.012. doi:10.1016/j.dci.2007.03.012,

  35. 35.

    Molina-Cruz A, DeJong RJ, Charles B, Gupta L, Kumar S, Jaramillo-Gutierrez G, Barillas-Mury C: Reactive oxygen species modulate anopheles gambiae immunity against bacteria and plasmodium. J Biol Chem. 2008, 283 (6): 3217-3223. 10.1074/jbc.M705873200. doi:10.1074/jbc.M705873200,

  36. 36.

    Felix R, Silveira H: The role of anopheles gambiae p450 cytochrome in insecticide resistance and infection. Insecticides - Pest Engineering. 2012, Rijeka, Crotia: InTech,

  37. 37.

    Vlachou D, Schlegelmilch T, Christophides GK, Kafatos FC: Functional genomic analysis of midgut epithelial responses in anopheles during plasmodium invasion. Curr Biol: CB. 2005, 15 (13): 1185-1195. 10.1016/j.cub.2005.06.044. doi:10.1016/j.cub.2005.06.044,

  38. 38.

    Baton LA, Robertson A, Warr E, Strand MR, Dimopoulos G: Genome-wide transcriptomic profiling of anopheles gambiae hemocytes reveals pathogen-specific signatures upon bacterial challenge and plasmodium berghei infection. BMC Genomics. 2009, 10 (1): 257-10.1186/1471-2164-10-257. doi:10.1186/1471-2164-10-257. Accessed 2013-04-19,

  39. 39.

    Dijkers PF, O’Farrell PH: Drosophila calcineurin promotes induction of innate immune responses. Curr Biol. 2007, 17 (23): 2087-2093. 10.1016/j.cub.2007.11.001. doi:10.1016/j.cub.2007.11.001. Accessed 2013-04-16,

  40. 40.

    Wähe A, Kasmapour B, Schmaderer C, Liebl D, Sandhoff K, Nykjaer A, Griffiths G, Gutierrez MG: Golgi-to-phagosome transport of acid sphingomyelinase and prosaposin is mediated by sortilin. J Cell Sci. 2010, 123 (14): 2502-2511. 10.1242/jcs.067686. doi:10.1242/jcs.067686. Accessed 2013-04-16,

  41. 41.

    Franc NC, Heitzler P, Alan BR, Ezekowitz RAB, White K: Requirement for croquemort in phagocytosis of apoptotic cells in drosophila. Science. 1999, 284 (5422): 1991-1994. 10.1126/science.284.5422.1991. doi:10.1126/science.284.5422.1991. Accessed 2013-04-19,

  42. 42.

    Gregorio ED, Spellman PT, Rubin GM, Lemaitre B: Genome-wide analysis of the drosophila immune response by using oligonucleotide microarrays. Proc Natl Acad Sci U S A. 2001, 98 (22): 12590-12595. 10.1073/pnas.221458698. doi:10.1073/pnas.221458698. Accessed 2013-11-12,

  43. 43.

    Scharlaken B, De Graaf DC, Memmi S, Devreese B, Van Beeumen J, Jacobs FJ: Differential protein expression in the honey bee head after a bacterial challenge. Arch Insect Biochem Physiol. 2007, 65 (4): 223-237. 10.1002/arch.20179. doi:10.1002/arch.20179. Accessed 2013-11-13,

  44. 44.

    Buttstedt A, Moritz RFA, Erler S: Origin and function of the major royal jelly proteins of the honeybee (Apis mellifera) as members of the yellow gene family. Bio Rev Camb Philos Soc. 2013, 89 (2): 255-269. doi:10.1111/brv.12052. Accessed 2013-11-13,

  45. 45.

    Amdam GV, Simões ZLP, Hagen A, Norberg K, Schroder K, Mikkelsen O, Kirkwood TBL, Omholt SW: Hormonal control of the yolk precursor vitellogenin regulates immune function and longevity in honeybees. Exp Gerontol. 2004, 39 (5): 767-773. 10.1016/j.exger.2004.02.010. doi:10.1016/j.exger.2004.02.010,

  46. 46.

    Korner P, Schmid-Hempel P: In vivo dynamics of an immune response in the bumble bee bombus terrestris. J Invertebr Patho. 2004, 87: 59-66. 10.1016/j.jip.2004.07.004.

  47. 47.

    van Rij RP, Saleh M-C, Berry B, Foo C, Houk A, Antoniewski C, Andino R: The RNA silencing endonuclease argonaute 2 mediates specific antiviral immunity in drosophila melanogaster. Genes Dev. 2006, 20: 2985-2995. 10.1101/gad.1482006. doi:10.1101/gad.1482006,

  48. 48.

    Gerardo NM, Altincicek B, Anselme C, Atamian H, Barribeau SM, Duncan EJ, Evans JD, Gabaldón T, Ghanim M, Heddi A, Kaloshian I, Latorre A, Moya A, Nakabachi A, Parker BJ, Pérez-Brocal V, Pignatelli M, Ramsey JS, Spragg CJ, Tamames J, Tamarit D, Tamborindeguy C, Vincent-Monegat C, Vilcinskas A, Vos Md: Immunity and other defenses in pea aphids, acyrthosiphon pisum. Genome Biol. 2010, 11 (2): 21-10.1186/gb-2010-11-2-r21. doi:10.1186/gb-2010-11-2-r21. Accessed 2013-04-19,

  49. 49.

    Ferdig MT, Taft AS, Smartt CT, Lowenberger CA, Li J, Zhang J, Christensen BM: Aedes aegypti dopa decarboxylase: gene structure and regulation. Insect Mol Biol. 2000, 9 (3): 231-239. 10.1046/j.1365-2583.2000.00187.x.

  50. 50.

    Ramet M, Manfruelli P, Pearson A, Mathey-Prevot B, Ezekowitz RAB: Functional genomic analysis of phagocytosis and identification of a drosophila receptor for e-coli. Nature. 2002, 416: 644-648. 10.1038/nature735.

  51. 51.

    Tahar R, Boudin C, Thiery I, Bourgouin C: Immune response of anopheles gambiae to the early sporogonic stages of the human malaria parasite plasmodium falciparum. EMBO J. 2002, 21 (24): 6673-6680. 10.1093/emboj/cdf664. doi:10.1093/emboj/cdf664. Accessed 2013-04-23,

  52. 52.

    Cirimotich CM, Dong Y, Garver LS, Sim S, Dimopoulos G: Mosquito immune defenses against plasmodium infection. Dev Comparative Immunol. 2010, 34 (4): 387-395. 10.1016/j.dci.2009.12.005. doi:10.1016/j.dci.2009.12.005. Accessed 2014-02-10,

  53. 53.

    Ligoxygakis P, Bulet P, Reichhart J-M: Critical evaluation of the role of the toll-like receptor 18-wheeler in the host defense of drosophila. EMBO Rep. 2002, 3 (7): 666-673. 10.1093/embo-reports/kvf130. doi:10.1093/embo-reports/kvf130. Accessed 2013-04-24,

  54. 54.

    Bond D, Foley E: A quantitative RNAi screen for JNK modifiers identifies pvr as a novel regulator of drosophila immune signaling. PLoS Pathog. 2009, 5 (11): 1000655-10.1371/journal.ppat.1000655. doi:10.1371/journal.ppat.1000655,

  55. 55.

    Corby-Harris V, Drexler A, Watkins de Jong L, Antonova Y, Pakpour N, Ziegler R, Ramberg F, Lewis EE, Brown JM, Luckhart S, Riehle MA: Activation of akt signaling reduces the prevalence and intensity of malaria parasite infection and lifespan in anopheles stephensi mosquitoes. PLoS Pathog. 2010, 6 (7): 1001003-10.1371/journal.ppat.1001003. doi:10.1371/journal.ppat.1001003. Accessed 2014-02-10,

  56. 56.

    Kocks C, Cho JH, Nehme N, Ulvila J, Pearson AM, Meister M, Strom C, Conto SL, Hetru C, Stuart LM, Stehle T, Hoffmann JA, Reichhart J-M, Ferrandon D, Rämet M, Ezekowitz RAB: Eater, a transmembrane protein mediating phagocytosis of bacterial pathogens in drosophila. Cell. 2005, 123 (2): 335-346. 10.1016/j.cell.2005.08.034. doi:10.1016/j.cell.2005.08.034. Accessed 2013-04-24,

  57. 57.

    Arakane Y, Muthukrishnan S, Beeman RW, Kanost MR, Kramer KJ: Laccase 2 is the phenoloxidase gene required for beetle cuticle tanning. Proc Natl Acad Sci U S A of the U S A. 2005, 102 (32): 11337-11342. 10.1073/pnas.0504982102. doi:10.1073/pnas.0504982102. Accessed 2014-02-13,

  58. 58.

    Blumberg BJ, Trop S, Das S, Dimopoulos G: Bacteria- and IMD pathway-independent immune defenses against plasmodium falciparum in anopheles gambiae. PLoS ONE. 2013, 8 (9): doi:10.1371/journal.pone.0072130. Accessed 2014-02-10,

  59. 59.

    Vierstraete E, Verleyen P, Baggerman G, D’Hertog W, Arckens L, Loof AD, Schoofs L, Bergh GVd: A proteomic approach for the analysis of instantly released wound and immune proteins in drosophila melanogaster hemolymph. Proc Natl Acad Sci U S A. 2004, 101 (2): 470-475. 10.1073/pnas.0304567101. doi:10.1073/pnas.0304567101. Accessed 2013-04-23,

  60. 60.

    Mukherjee T, Kim WS, Mandal L, Banerjee U: Interaction between notch and hif-alpha in development and survival of drosophila blood cells. Science (New York, N.Y.). 2011, 332 (6034): 1210-1213. 10.1126/science.1199643. doi:10.1126/science.1199643,

  61. 61.

    Watson FL, Püttmann-Holgado R, Thomas F, Lamar DL, Hughes M, Kondo M, Rebel VI, Schmucker D: Extensive diversity of ig-superfamily proteins in the immune system of insects. Science. 2005, 309 (5742): 1874-1878. 10.1126/science.1116887. doi:10.1126/science.1116887. Accessed 2013-04-19,

  62. 62.

    Dong YM, Taylor HE, Dimopoulos G: AgDscam, a hypervariable immunoglobulin domain-containing receptor of the anopheles gambiae innate immune system. Plos Biol. 2006, 4: 1137-1146.

  63. 63.

    Dong Y, Cirimotich CM, Pike A, Chandra R, Dimopoulos G: Anopheles NF-kappaB-regulated splicing factors direct pathogen-specific repertoires of the hypervariable pattern recognition receptor AgDscam. Cell Host Microbe. 2012, 12 (4): 521-530. 10.1016/j.chom.2012.09.004. doi:10.1016/j.chom.2012.09.004,

  64. 64.

    Mayans O, Benian GM, Simkovic F, Rigden DJ: Mechanistic and functional diversity in the mechanosensory kinases of the titin-like family. Biochem Soc Trans. 2013, 41 (4): 1066-1071. 10.1042/BST20130085. doi:10.1042/BST20130085,

  65. 65.

    Garver LS, Xi Z, Dimopoulos G: Immunoglobulin superfamily members play an important role in the mosquito immune system. Dev Comp Immunol. 2008, 32 (5): 519-531. 10.1016/j.dci.2007.09.007. doi:10.1016/j.dci.2007.09.007,

  66. 66.

    Castro DP, Moraes CS, Gonzalez MS, Ratcliffe NA, Azambuja P, Garcia ES: Trypanosoma cruzi immune response modulation decreases microbiota in rhodnius prolixus gut and is crucial for parasite survival and development. PLoS ONE. 2012, 7 (5): 36591-10.1371/journal.pone.0036591. doi:10.1371/journal.pone.0036591. Accessed 2013-07-24,

  67. 67.

    Oberle M, Balmer O, Brun R, Roditi I: Bottlenecks and the maintenance of minor genotypes during the life cycle of trypanosoma brucei. PLoS Pathog. 2010, 6 (7): 1001023-10.1371/journal.ppat.1001023. doi:10.1371/journal.ppat.1001023. Accessed 2014-04-28,

  68. 68.

    Beckerman AP, de Roij J, Dennis SR, Little TJ: A shared mechanism of defense against predators and parasites: chitin regulation and its implications for life-history theory. Ecol Evol. 2013, 3 (15): 5119-5126. 10.1002/ece3.766. doi:10.1002/ece3.766,

  69. 69.

    Weiss BL, Wang J, Maltz MA, Wu Y, Aksoy S: Trypanosome infection establishment in the tsetse fly gut is influenced by microbiome-regulated host immune barriers. PLoS Pathog. 2013, 9 (4): 1003318-10.1371/journal.ppat.1003318. doi:10.1371/journal.ppat.1003318. Accessed 2013-10-03,

  70. 70.

    Mallon EB, Loosli R, Schmid-Hempel P: Specific versus nonspecific immune defense in the bumblebee, bombus terrestris l. Evolution. 2003, 7: 1444-1447.

  71. 71.

    Moret Y, Schmid-Hempel P: Immune responses of bumblebee workers as a function of individual and colony age: senescence versus plastic adjustment of the immune function. Oikos. 2009, 118: 371-378. 10.1111/j.1600-0706.2008.17187.x.

  72. 72.

    Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, Robles M, Talón M, Dopazo J, Conesa A: High-throughput functional annotation and data mining with the blast2go suite. Nucleic Acids Res. 2008, 36 (10): 3420-3435. 10.1093/nar/gkn176. doi:10.1093/nar/gkn176. Accessed 2013-03-27,

  73. 73.

    Wu TD, Nacu S: Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics (Oxford, England). 2010, 26 (7): 873-881. 10.1093/bioinformatics/btq057. doi:10.1093/bioinformatics/btq057,

  74. 74.

    McCarthy DJ, Chen Y, Smyth GK: Differential expression analysis of multifactor RNA-seq experiments with respect to biological variation. Nucleic Acids Res. 2012, 40 (10): 4288-4297. 10.1093/nar/gks042. doi:10.1093/nar/gks042,

  75. 75.

    Team R-C: R: A language and environment for statistical computing. 2013, Vienna, Austria: R foundation for statistical computing

  76. 76.

    Supek F, Bošnjak M, Škunca N, Šmuc T: REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE. 2011, 6 (7): 21800-10.1371/journal.pone.0021800. doi:10.1371/journal.pone.0021800. Accessed 2013-04-16,

  77. 77.

    Brucker RM, Funkhouser LJ, Setia S, Pauly R, Bordenstein SR: Insect innate immunity database (IIID): An annotation tool for identifying immune genes in insect genomes. PLoS ONE. 2012, 7 (9): 45125-10.1371/journal.pone.0045125. doi:10.1371/journal.pone.0045125. Accessed 2013-09-23,

  78. 78.

    Trapnell C, Pachter L, Salzberg SL: TopHat: discovering splice junctions with RNA-seq. Bioinformatics. 2009, 25 (9): 1105-1111. 10.1093/bioinformatics/btp120. doi:10.1093/bioinformatics/btp120. Accessed 2014-05-27,

  79. 79.

    Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L: Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012, 7 (3): 562-578. 10.1038/nprot.2012.016. doi:10.1038/nprot.2012.016. Accessed 2014-02-26,

  80. 80.

    Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup: The sequence alignment/map format and SAMtools. Bioinformatics (Oxford, England). 2009, 25 (16): 2078-2079. 10.1093/bioinformatics/btp352. doi:10.1093/bioinformatics/btp352,

  81. 81.

    Anders S, Pyl PT, Huber W: HTSeq - a python framework to work with high-throughput sequencing data. bioRxiv. 2014, doi:10.1101/002824. Accessed 2014-05-27,

  82. 82.

    Anders S, Reyes A, Huber W: Detecting differential usage of exons from RNA-seq data. Genome Res. 2012, 22 (10): 2008-2017. 10.1101/gr.133744.111. doi:10.1101/gr.133744.111,

  83. 83.

    Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local alignment search tool. J Mol Biol. 1990, 215 (3): 403-410. 10.1016/S0022-2836(05)80360-2. doi:10.1006/jmbi.1990.9999,

Download references


Thanks to Simon Anders and Alejandro Reyes for help with the DEXSeq model Thanks to Paul Schmid-Hempel for discussions. CR was funded by a BBSRC studentship. This work was partially funded through a NERC NBAF pilot grant (NBAF606) to EBM. We thank the Bumblebee Genome Consortium ( for providing genomic resources that were used for this study.

Author information

Correspondence to Eamonn B Mallon.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

CER carried out the infections and extracted RNA. JDLG prepared samples for library production. SA analysed the data and wrote the manuscript. SMB analysed the data. DT wrote the manuscript. EBM analysed the data and wrote the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material


Additional file 1: Table S1: The sequences, statistics and annotations for all differentially expressed genes in each of the three differential expression analyses are available from (XLSX 632 KB)

Additional file 2: Table S2: The complete output of the DEXSeq analysis is available from (ZIP 44 KB)


Additional file 3: Table S3: 615 loci displayed alternative splicing depending on the the interaction between the host genotype and the parasite genotype. The sequences, statistics and annotations for these loci is available in the supplementary data from (XLSX 4 MB)

Authors’ original submitted files for images

Rights and permissions

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark


  • Genotype-genotype
  • Peritrophic membrane
  • Social insects