Skip to main content

In the battle of the disease: a transcriptomic analysis of European foulbrood-diseased larvae of the Western honey bee (Apis mellifera)



European foulbrood is a significant bacterial brood disease of Apis sp. and can cause severe and devastating damages in beekeeping operations. Nevertheless, the epidemiology of its causative agent Melissococcus plutonius has been begun to uncover but the underlying mechanisms of infection and cause of disease still is not well understood. Here, we sought to provide insight into the infection mechanism of EFB employing RNAseq in in vitro reared Apis mellifera larvae of two developmental stages to trace transcriptional changes in the course of the disease, including Paenibacillus alvei secondary infected individuals.


In consideration of the progressing development of the larva, we show that infected individuals incur a shift in metabolic and structural protein-encoding genes, which are involved in metabolism of crucial compounds including all branches of macronutrient metabolism, transport protein genes and most strikingly chitin and cuticle associated genes. These changes underpin the frequently observed developmental retardation in EFB disease. Further, sets of expressed genes markedly differ in different stages of infection with almost no overlap. In an earlier stage of infection, a group of regulators of the melanization response cascade and complement component-like genes, predominantly C-type lectin genes, are up-regulated while a differential expression of immune effector genes is completely missing. In contrast, late-stage infected larvae up-regulated the expression of antimicrobial peptides, lysozymes and prominent bacteria-binding haemocyte receptor genes compared to controls. While we clearly show a significant effect of infection on expressed genes, these changes may partly result from a shift in expression timing due to developmental alterations of infection. A secondary infection with P. alvei elicits a specific response with most of the M. plutonius associated differential immune effector gene expression missing and several immune pathway genes even down-regulated.


We conclude that with progressing infection diseased individuals undergo a systemic response with a change of metabolism and their activated immune defence repertoire. Moreover, larvae are capable of adjusting their response to a secondary invasion in late stage infections.

Peer Review reports


Uncovering disease mechanisms is crucial to get a deeper insight into infection dynamics and ecological impact of diseases on host population development and host-pathogen interactions. Most important, acquired knowledge contributes to a broader understanding of the implications and conditions in epidemiological scenarios, and enables the development of strategies for disease control on different levels of interventions and appropriate measures to counteract or even prevent disease outbreaks. European foulbrood (EFB) is a widely distributed and significant honey bee brood disease but despite more than one century of research there is still insufficient understanding of this disease on molecular level [1, 2]. Recently, large regional outbreaks of EFB have been reported from e.g. England and Wales or Switzerland and many cases across numerous European countries were recorded and confirmed, and continue to persistently re-emerge in affected regions [2]. Other global areas like the United States also observe significant increase of EFB cases (pers. comm. Marla Spivak). Changes in beekeeping practises, agricultural intensification [3] and expanding global trade with worldwide accelerated exchange of stock [4] (and not at least climate change), raise presentiments of an upsurge of new epidemics and endemic outbreaks. This causes a pressing need to shed light on the basic mechanisms, the course and the most significant factors contributing to severe infections and disease.

Besides an urgency of developing efficient programs and treatment methods to prevent and counteract new outbreaks, several features of the Gram-positive bacterium Melissococcus plutonius, the causative agent of EFB, make it an interesting and promising model pathogen to illuminate general pathogenic mechanisms and counteracting host immune responses ranging from the basic cellular to a behavioural colony level [5, 6]. Its features are the wide range of virulence diversity of different strains [7, 8] and regional enzootic states [9], its resilience to desiccation, high sugar concentrations, and spore-independent viability for months as well as a brood specialization [2].

M. plutonius typically infects honey bee larvae in the first days after hatching via ingestion of contaminated food. Then, the larval susceptibility to infection and development of severe disease decreases with larval age until there is no pathogenic effect in the final larval stage before pupation [1, 2]. This high susceptibility of young larvae is also observed in another severe brood diseases, i.e. American Foulbrood (AFB) [2], caused by Paenibacillus larvae, suggesting a similar, host-dependent cause of limited infectivity (or a pathogen barrier) in the late larval stages for EFB and AFB, respectively. Although, different mechanisms of resistance and implicated factors have been discussed in former studies, the most critical factor seems to be the pathogen multiplication time until pupation with a potential clearance of the pathogen and a lag of bacterial toxin release [8, 10]. Former studies reported presumably EFB-resistant (or rather tolerant) larvae found in naturally as well as experimentally infected colonies, an observation we also regularly made in our in vitro studies [6,7,8, 11,12,13]. These larvae survive pathogen infestation, succeed to pupate (after defecation) and finally emerge as adult workers. A fraction of these workers is seemingly underweight but some of them are morphologically comparable to non-infected individuals. Infections usually depend not solely on the inoculation doses of the infective agent but on many factors like timing of infection, the genotype of the pathogen, the actual susceptibility of the host represented by its genotype and constitution, and also on stochastic processes [14, 15].

Different hypotheses for the observed defence in larvae have been considered, encompassing any layer of immune defence including e.g. the variation in physical barrier properties of the gut, immune system activity comprising immune effector gene expression differences in larvae [7, 13, 16]. Further, food composition and worker added substances have been brought forward as potential colony level defences conveying individual resistance of the brood [17]. Aside from that, a mechanism was suggested where an infestation of larvae leads to a host-pathogen competition for nutrients [11, 12]. In this scenario, the multiplication of the bacterium inside the food mass leads to potential malnutrition, developmental retardation and ultimately starvation and death of the infested larva like hypothesised and substantiated by histological examinations [18,19,20]. Starvation of infected larvae also indicates early stages of P. larvae infections. This was implied by a proteomic study in which multiplying vegetative cells of P. larvae metabolize glucose and fructose, the main carbohydrates of food jelly, without detectable tissue damage [21]. Although, M. plutonius seems to adhere to the peritrophic matrix (PM) it is clearly contained inside the gut lumen and does not breach the gut epithelium barrier until the larva dies and decay starts [20, 22]. This characteristic course is distinctly different from the disease mechanism of AFB. In AFB the causative agent P. larvae actively attacks and breaches the PM structure and the ventricular epithelium, and ultimately invades the haemocoel of the infected larvae [23]. This elicits the typical immune response observed in AFB affected larvae while a response in M. plutonius infection was not demonstrated so far [5, 16].

Here, we present an experimental in vitro approach to get a general overview of the expressional response of honey bee larvae artificially infected with M. plutonius and a common secondary invader, Paenibacillus alvei. This enables an onset of a deeper understanding of the EFB disease infection mechanism and the immune response on a physiological level.


Bacterial infections may elicit complex and profound systemic responses with far reaching consequences for the attacked organisms. In honey bee larvae, infection with a pathogenic strain of M. plutonius often leads to a retardation of growth and therefore interfering with larval metabolism and development with a fatal outcome in severe infections [1].

A principal component analysis (PCA) revealed that all larvae samples from the same treatment groups clustered together after batch effect correction. The data set was strongly separated by larval age (day 4 = 3 days post infection (dpi) vs. day 7 = 6 dpi, 42.53% of total variance) in one dimension, rather than by infection (12.66% of total variance) in the young larvae, compared to the second sampling date where older larvae were clearly separated in distinct clusters (Fig. 1). These results illustrate that expressional developmental changes over time elicit a stronger effect than treatment, i.e. infection. Interestingly, we also find separation of clusters in our correlation of larval weight and M. plutonius plasmid copy number (Fig. S 1, Additional file 1). Moreover, comparing DEG sets of contrasts of day 3 (Cont vs. Inf on 3 dpi) and day 6 post infection (Cont vs. Inf; Cont vs. secInf on 6 dpi), a strong difference of expressional patterns was evident with almost no overlap of DEGs between both sampling time points (FDR < 0.05, − 0.585 < log2FC > 0.585; Fig. 2).

Fig. 1
figure 1

PCA with filtered data set (batch effects removed by RUVr method from the RUVSeq R package). Data sets clearly cluster in respective groups and are markedly separated in larval age groups (3 dpi vs. 6 dpi), while a strong pronounced separation of infected and control samples without an overlap was only observed for 6 dpi. Numbers refer to replicate of group. For more details see Table S 1 (Additional file 3) and Fig. S 8 (Additional file 1) where colors correspond to the same groups, respectively. M.p. – M. plutonius, P.a. – P. alvei

Fig. 2
figure 2

Overlap of single group contrasts of control and infected larvae. a Down-regulated DEGs and b Up-regulated DEGs. Expressional differences with a high fold change (FC > > 2). (Cont: control, Inf: M. plutonius-infected, secInf: M. plutonius infected + P.alvei)

In a functional annotation clustering (FAC) approach, we further focused on the developmental transition in mRNA abundances in the treatment groups. We first addressed and illustrated the transcriptional change between two time points in the larval development (3 dpi and 6 dpi) and compared the respective differentially expressed gene (DEG) set contrast of 3 dpi vs. 6 dpi of the controls and the contrast of 3 dpi vs. 6 dpi of the M. plutonius infected. Then we focused on mRNA abundance of candidate genes of transport/storage proteins and innate immune system genes, which were shown to be a part of a general response in other infections [5, 21]. Finally, we analysed the effect of a secondary infection with P. alvei.

Correlation of weight and expressed genes

Overall, transcriptional activity was increased in lighter and hence younger larvae where mRNA levels of 1110 genes were positively and 2310 negatively significantly correlated with larval weight (Fig. S 2, Additional file 1; Table S 1, Additional file 3). For the positively weight correlated genes, highest ranking terms in the FAC analysis were associated with cellular transport, endocytosis (GO:0035091, GO:0016192, GO:0006886, GO:0035091, ame04144) and lysosome (ame04142) as well as fatty acid degradation (ame00071). In contrast, the negatively weight correlated genes and therefore stronger expressed in young and lighter larvae were mostly involved in splicing, translation and in general RNA turnover associated terms, most strikingly represented by a large set of ribosome protein-encoding genes (ame03010) with consistent correlation coefficients. The other large group of expressed genes is involved in oxidative phosphorylation (ame00190) (for additional details see Table S 1, Additional file 3).

Functional annotation clustering in temporal expressional shift

As a general overview of transcriptional change in controls and infected larvae, we performed an Apis mellifera specific FAC analysis to identify enriched terms and retrieve functional information of enriched genes from the DEG sets (Table S 2, Additional file 4). In total, we identified 3070 DEGs (p < 0.05, FDR-corrected: − 1 < log2FC > 1) which were significantly regulated. 2090 of the 2255 DEGs in the control group and 1712 of 1848 DEGs in the infected group (Table 1, Fig. S 3, Additional file 1), with A. mellifera specific accession in the DAVID database and a differential expression over time (3 dpi vs. 6 dpi, respectively; 1251 DEGs overlapping between both contrasts), were included in our FAC analysis (Table S 2, Additional file 4). FAC yielded 9 annotation clusters containing significant terms (87 in total) for the change of mRNA levels of the controls and 15 annotation clusters (73 in total) for the infected group (p < 0.1; Table S 2, Additional file 4). For both controls and infected the two highest-ranking clusters were associated with cuticle proteins, extracellular localisation and chitin binding, illustrating the importance in larval development and the body and gut lining (Table 2; Fig. 3). This was also substantiated by the other categories where the highest-ranking significant terms from Molecular Function (GO:0042302, GO:0008061) and Biological Processes (GO:0006030, GO:0006633) were associated with chitin binding and metabolism. Chitin-binding proteins (CBP) like perithrophins and cuticular proteins but also mucins are a constitutive part of these extracellular structures. In our study, several CBPs undergo significant differential regulation in infected larvae compared to controls (Table S 3, Additional file 5). Corresponding, in the overlap of gene ontology (GO) terms of controls and M. plutonius infected the significant terms for “Cellular localisation” were identified as the extracellular region and the membrane (GO:0005576, GO:0016021; Table 2).

Table 1 Group contrasts of significantly differentially expressed genes and transcripts
Table 2 Functional annotation clustering summary of temporal expressional change
Fig. 3
figure 3

Enrichment network of the DAVID functional annotation analysis. Terms are represented by circles and overlap of terms is represented by node connection resulting in clusters of terms with overlapping functions. Circle size is related to enriched gene set size and edge width is positively associated with overlap of gene sets between the respective terms. For detailed enrichment network see Fig. S 11 (Additional file 2)

Metabolism related transcriptional changes (3 dpi vs. 6 dpi), particularly represented by Kyoto Encyclopedia of Genes and Genomes (KEGG) terms, significantly enriched in both control and infected groups, involved all major branches like amino acid, carbohydrate and lipid metabolism. Although, most pathways overlapped between controls and infected individuals, there were some terms enriched uniquely and/or expressed exclusively for one or the other group (Table 2). In infected larvae several pathways were differentially regulated on expressional level. In general, in the infected group (3 dpi vs. 6 dpi) just 109 genes (154 in controls) of the parent term “Metabolic pathways” (ame01100) were significantly regulated over time with just 84 genes overlapping in infected and controls and 25 genes showing significant expression in the infected group exclusively (Table 3). This points to a disruption or temporal shift of metabolic pathway regulation in infected larvae. The differentially enriched term “Integral component of membrane” (GO:0016021) is connected to metabolism by the pronounced role of membrane located proteins in lipid metabolism (GO:0006633 = fatty acid biosynthetic process) but also crucial substrate transporter proteins (Table S 2, Additional file 4).

Table 3 Significantly enriched KEGG and InterPro terms of DAVID functional annotation in temporal expressional change

In the course of infection (Inf 3 dpi vs. Inf 6 dpi) significantly expressed gene sets of the carbon metabolism (ame01200) as well as the overlapping pentose phosphate pathway (ame00030) and glycolysis/gluconeogenesis (ame00010) are strongly reduced compared to controls (Cont 3 dpi vs. Cont 6 dpi). In contrast, the expressional pattern of the pyruvate pathway (ame00620) gene set composition seems to be sustained for the most part. All these pathways are crucial in storage and conversion of energy reserves for development and metamorphosis in insect larvae. Moreover, there are two functionally related InterPro terms of transporter proteins with a differential enrichment, the “General substrate transporters” (IPR005828) and the “Major facilitator superfamily” (IPR020846) (Table 3). These terms comprise two overlapping sets of membrane transport proteins of various molecules but most importantly including carbohydrate and other organic molecule transporters. Specifically, the change in the trehalose and glucose transporter expression points to a shift in carbohydrates conversion, storage and distribution in infected larvae over time (for additional details see Table S 2, Additional file 4).

Temporal expression changes of several genes from different classes of protein families (InterPro category) were significantly enriched in the infected larvae particularly Trypsin-like cysteine/serine like peptidase (IPR009003), CRAL/TRIO domain (IPR001251), heat shock protein 20 (Hsp20) domain (IPR002068), DUF1676 (IPR012464) and Zona pellucida (ZP) domain (IPR001507) containing protein genes (Fig. 3, Table 3; for additional details Table S 2, Additional file 4).

Specific proteins involved in larval development

Insect larvae possess specific storage and transport proteins. In honey bees, proteins like hexamerines and a Very high density lipoprotein (VHDL) are strongly expressed in larvae and are accumulated in the haemolymph and the fat body preceding metamorphosis, respectively [24, 25]. Other functional proteins, like fibroins are massively synthesised in late-stage larvae prior to the spinning stage and are indicative for proper developmental progress [26]. We found a significant expressional change of the corresponding protein-encoding genes when comparing young 3 dpi and old 6 dpi larvae in both, infected and control groups, respectively (Fig. 4A). While expression of genes encoding storage and transport proteins was mostly unchanged in infection with M. plutonius as well as in secondary infected larvae (M. plutonius + P. alvei treatment), silk fibroin mRNA levels were massively down-regulated with a factor of 24 to more than 40 in infected larvae. Notably, secondary infection (M. plutonius + P. alvei) led to a similar reduced silk fibroin expression in old larvae while storage protein mRNA stayed unchanged and even dropped for vhdl and therefore was failed to be regulated (Fig. 4A).

Fig. 4
figure 4

Candidate gene expression of a larval specific and storage/transport protein genes and b canonical immune effector and signalling genes in single group comparisons. (Dark green: down-regulation, FC > 1; light green: down-regulation, FC < 1; purple: up-regulation, FC > 1; pink: up-regulation, FC < 1; light grey: not differentially regulated; dark grey: not detected). In case of a borderline significance the p-value is shown for the respective comparison. Expressional differences with a high fold change (FC > > 2) but not significant values are marked with N.S. For more details see Table S 6 (Additional file 8)

In our FAC analysis, we found several genes encoding groups of proteins, which were solely enriched in temporal expressional change of infected larvae (Table 4). DUF1676 is characteristic for the osiris gene family which were found to be coincidentally expressed with cuticle deposition in development but are also suspected to be involved in other functions like immunity [27]. We found a number of osiris genes to be significantly regulated over time in the control group (3 dpi vs 6 dpi) and even more pronounced in the infected larvae of up to a 1000-fold change. While at 6 dpi there were no differences between control and infected group, the 3 dpi contrast showed several osiris genes to be significantly up-regulated in the infected individuals (Table S 4, Additional file 6). The Zona pellucida (ZP) domain is associated with diverse functions but primarily is present in membrane-bound proteins with filaments and/or matrix related context [28]. In infected larvae, genes of ZP domain containing proteins were mostly down-regulated over time compared to the control group. All genes belonging to the stress responsive group of Hsp20 chaperones, including several l(2)efl genes, were up-regulated in older infected individuals (Table 4).

Table 4 Differentially enriched InterPro terms of functional groups from FAC analysis of significant DEGs

Innate immune response and immune pathway signalling

Toll and Imd/JNK signalling

Honey bee larvae express a strong immune response upon a challenge with pathogens from different types of organisms [5]. In our experimental set up, we therefore expected an activation of the Toll and Imd pathways, which are responsive to peptidoglycan containing Gram-negative and Gram-positive bacteria, like M. plutonius. In regard of the mRNA abundance of Toll and Imd signalling pathway related genes of the control groups (Cont 3 dpi vs. Cont 6 dpi), we observed a strong up-regulation of bacteria specific pattern recognition receptors (PRR) for both Toll and Imd pathway over time while a fungi specific PRR was significantly down-regulated (b-gluc1, FC = 2.71; Fig. 4B). Further, a decrease of the Imd inhibitor pirk and the dual oxidase expression was observed (Fig. 4B). While pirk is significantly down-regulated in control larvae over time (FC = 4.0, p = 0.0053), a temporal change in the infected groups was not recorded. (Inf 3 dpi vs. Inf 6 dpi; FC = 0.01, p = 0.99; Table S 5, Additional file 7).

Innate immunity

We detected only two differentially expressed canonical immunity genes in samples of day 3 post infection, phenol oxidase-activating factor 2 and proclotting enzyme SP7, serine proteases that were shown to be involved in prophenoloxidase (PPO) activation [29]. In addition, several complement component associated protein genes like C-type lectins (CTL1, CTL4, CTL6, CTL12) and another complement-like protein (LOC725964) were two- to four-fold up-regulated in infected individuals (Fig. 4B; Table S 6, Additional file 8).

In our experimental set up the onset of immune effector expression was observed with progressing infection in the 6 dpi larvae. We detected increased mRNA levels of several antimicrobial peptides (AMPs) like hymenoptaecin (7-fold) and abaecin (more than 2-fold) and lysozymes (LOC113218576, LOC724899; 3.5-fold and 2.5-fold) in older infected larvae (Fig. 4B). Defensins comprise a prominent and widely conserved AMP family and usually are expected to be expressed upon an immune challenge. In our data set defensin 1 and 2 where not represented or just inconsistently and in negligible read count numbers. The expressional differences of regulated AMPs (including apidaecin) in treatment groups are exceeded by a magnitude when compared to the differential expression of these AMPs over time (controls as well as infected). Here, in both controls and infected we found a strong increase of mRNA levels up to more than 6000 times (abaecin) higher than in 3 dpi larvae.

However, the most striking transcriptional change in 6 dpi infected larvae (both just M. plutonius and secondarily infected) was an increase in expression of cellular defence associated genes eater-like (LOC724421) and a scavenger receptor class B member 1 (LOC413408; also peste homologue) compared to 6 dpi controls (Fig. 4B). Particularly, eater-like is neither regulated in the controls nor infected larvae over time (3 vs 6 dpi) and therefore seems to be a specific expressional response to the induced infection.

Feeding behaviour and hormone-associated genes

In M. plutonius infections we usually observe diseased larvae with a strong reduction in body weight and cessation of feeding [6, 7]. Hence, we were interested in the expression of candidate genes related to feeding behaviour. Comparing 6 dpi controls and 6 dpi infected larvae we found differences in mRNA levels of genes, which previously have been associated with feeding behaviour and nutritional stress [30, 31]. Particularly, neuropeptide F (npf) was significantly higher expressed in 6 dpi infected larvae (3.2 times, FDR = 0.023) and even more pronounced in larvae secondary infected with P. alvei (5.3 times, FDR < 0.001) compared to controls. These results are substantiated by a differential expression of the serine protease scarface (scaf) gene, which was included in the group of IPR001254: Peptidase S1 term of our FAC analysis. Scaf was down-regulated in the 6 dpi M. plutonius infected group by a factor of 1.8 compared to 6 dpi controls (FDR < 0.001). In D. melanogaster Scaf acts as insulin signalling independent neuromodulator in sugar deprived states and is expressed in scaf-specific specialized neurons [31]. These data are further supported by our FAC of a transcriptional change in significantly enriched genes categorised in GO terms related to metabolism and transport of carbohydrates (Table S 2, Additional file 4). Interestingly, 6 dpi M. plutonius-infected larvae showed an altered expression of genes involved in hormone synthesis pathways (Table S 7, Additional file 9; Fig. S 4, Additional file 1). One of the important groups here are heme-containing enzyme genes. Several of these genes, which often possess oxidoreductase activity and may also be involved in detoxification of xenobiotics, were differentially expressed between infected and controls over time (3 dpi vs 6 dpi; Table S 8, Additional file 10). Specifically, altered expression of two farnesol dehydrogenases genes (LOC412458, LOC725489), enzymes involved in the juvenile hormone (JH) synthesis pathway, mediating the oxidation of farnesol to farnesal [32], may have implications for timing of larval development, successful metamorphosis and ultimately survival. Another important regulator of larval development is 20-Hydroxyecdysone (20-HE). A perturbation of the synthesis of active 20-HE through inhibition of the hydroxylation of ecdysone is indicated by down-regulated cytochrome P450 314A1 (also called shade) gene expression by a factor of two (Table S 7, Additional file 9; Fig. S 4, Additional file 1).

Secondary infection

We compared transcriptional levels of the DEG sets of solely M. plutonius infected 6 dpi (Cont vs Inf) with M. plutoinus and secondarily P. alvei infected (Cont vs secInf) 6 dpi larvae to determine the effect of secondary infection on gene expression patterns. We did not find significantly increased mortality in secondarily infected larvae in the recent study (compare [7, 8]) but, nevertheless, a marked difference of expressed genes of Cont vs Inf compared to Cont vs secInf individuals. Most DEGs from the Cont vs Inf contrast are also represented by comparable fold changes in the Cont vs secInf contrast, confirming the consistency of the DEG pattern of M. plutonius infection (Figs. 2 and 4). FAC analysis with a less restrictive DEG threshold (− 0.585 < log2FC > 0.585; p < 0.05; separated up- and down-regulated genes; for additional details see Table S 12, Additional file 14) yielded significant enrichment of 6 terms in the GO and KEGG categories for up-regulated DEGs in Cont vs Inf (348 DEGs) compared to 13 terms in Con vs secInf (696 DEGs). In secondarily infected individuals, most of these terms are functionally associated with metabolism, specifically a large set of KEGG terms related to carbohydrate metabolism and other terms with a function in release of stored energy like the pyruvate metabolism or the citric acid cycle. In contrast, in the solely M. plutonius infected group, the set of up-regulated terms is limited to galactose and glucose metabolism. Here, we found rather a more expanded set of down-regulated significant terms with function in transport, chitin binding and oxidoreductase activity (GO:0005215, GO:0008061, GO:0016614). As mentioned above, the response was for the most part comparable in the secondarily infected group but several DEGs were regulated differently (Fig. 4; Table S 6, Additional file 8). The lipoprotein genes of vitellogenin (vg) and apolipophorin-III-like protein (A4), and to a lesser amount the apolipophorin gene LOC408961, were significantly up-regulated in secInf compared to controls (fold change: vg = 3.4, A4 = 3.2, LOC408961 = 1.6). Surprisingly, in secondarily infected larvae mRNA levels of one lysozyme gene (LOC724899), abaecin and hymenoptaecin were not different from controls. Furthermore, there was a remarkable down-regulation of apidaecin (fold change: 7.1) in the secondarily infected group as well as a decreased expression level of two PRRs (fold change: PGRP-S3 = 2.0, LOC725832 = 2.8).

Regulatory non-coding RNAs

Non-coding RNAs (ncRNA) are involved in gene regulation in a variety of organisms and different types of regulatory RNAs were described in the last couple of decades [33]. Besides many other crucial functions, ncRNAs were shown to play an important role in innate immunity and are involved in pathogenic mechanisms of different pathogens [34]. While only one microRNA Mir9865 was significantly regulated in one group contrast, namely 3 dpi infected vs 6 dpi infected (p = 0.015, FC = 2.18, down-regulated), we found that regulation of several putative ncRNAs was strongly dependent on M. plutonius infection. In our expressed transcript set, we identified 1399 putative ncRNA specimen in both controls and infected larvae without being exclusive for one group. In total, we identified 8 ncRNAs which were significantly different in their in-group temporal change of abundance in controls versus infected (Table 5).

Table 5 Significant expressional change of non-coding RNAs putatively regulated in the course of M. plutonius-infection


We provide a transcriptome study on M. plutonius and P. alvei secondarily infected larvae of the Western honey bee Apis mellifera resembling European foulbrood disease with new insights on course of infection and its impact on larval development. First, we confirm previous findings that young larvae are highly active in nutrient conversion and protein synthesis with a high rate of oxidative phosphorylation, gene expression and translation [35]. Heavily enriched sets of genes involved in RNA-turnover and transport further confirm these results (Table S 1, Additional file 3). In older individuals, conversion and allocation of synthesised and stored reserves, including fatty acid degradation, are predominantly marked by processes related to inter- and intracellular trafficking, reflecting the high activity of the fat body cells and preparation to metamorphosis. In the next step of our analysis, we compared the temporal change of functional enriched gene sets of infected and control larvae employing functional annotation clustering.

Most notably, functional annotation clustering (FAC) analysis yielded significant enrichment of chitin metabolism and chitin biosynthesis-associated terms, ranking highest in the categories “Molecular Function” and “Biological Process”. These clusters contain a large group of peritrophines and cuticular protein genes as well as genes of other chitin binding proteins. Chitin turnover occupies a predominant role in larval development as chitin containing structures are repeatedly built up, broken down and recycled to enable larval growth and to constitute the body surface as well as the peritrophic matrix, an important component of the larval ventriculus [36]. Other structural proteins like collagens and laminins are important components of extracellular matrices (ECM), e.g. the lining of organs, and are involved in ECM receptor signalling and cell-to-cell communication [37]. Hence, perturbation of chitin and extracellular structure synthesis and metabolism will most likely affect the development of the larva.

A differential set of genes involved in several branches of metabolism indicate a change and shift in metabolic activity. Gene sets of several pathways presumably failed to be regulated in the same manner in the development (i.e. from 3 dpi to 6 dpi) of infected larvae as compared to control larvae. Here, carbon metabolism-associated pathways, which are involved in the metabolism of fat and glycogen seem to be strongly affected. This is also reflected in a change of associated membrane transport proteins expression. Larvae build up lipid and glycogen stores in fat body cells during the feeding stages and these are increasingly found in the larval trophocytes of the fat body [38]. In general, these resources are accumulated for energy supply and synthesis of structural and metabolic components, which are crucial in metamorphosis. An intruding pathogen elicits an immune system activation with an allocation of these resources [39]. Mounting an immune response is usually predicted to be costly and to trade off with other fundamental functions most strikingly seen under stress conditions like starvation and poor nutrient supply [15]. Whether the developmental retardation in EFB disease is a consequence of allocation of resources to the immune system or a direct result of competition for nutrients with the pathogen needs further investigations.

Two other groups of genes were represented in the FAC analysis and are associated with ECM, DUF1767 domain containing protein-encoding genes (predominantly Osiris family) and Zona pellucida (ZP) domain family genes. Although, their functions are not well understood, these results overlap with transcriptome studies of infections with another honey bee brood pathogen, namely P. larvae [40]. The most pronounced similarity is the enrichment and differential expression of Chitin-binding Protein (CBP) as well as Osiris family and ZP domain protein genes. This is not surprising as these infections are both located in the gut of the larva and most likely interfere with the larval development as already suspected by Cornman et al. [40]. Hence, this expressional pattern may be a part of a more general mechanism in infections of honey bee larvae. A shift or disturbance in larval development may also be derived from the observed lowered mRNA levels of the fibroins in M. plutonius infected and P. alvei secondary infected larvae in our study. Silk fibroin is synthesised prior to the spinning stage and used in the beginning of pupation [26]. Beyond that, the observed enrichment of a set of Hsp20 genes may be an indicator of a stress response in infected individuals in the course of the infection. Hsp20 proteins are a family of chaperones, which play a significant role in different physiological processes including the maintenance of cellular proteins under a variety of stress conditions [41].

Toll and Imd/JNK signalling

The observed decrease of Imd inhibitor pirk and the dual oxidase expression in untreated larvae is most likely attributable to an immune system activation in preparation to metamorphosis, in which pirk transcription decreases and therefore gives rise to immune effector expression like e.g. observed for the elevation of several AMP mRNA levels in the same group contrast. Coincidently, dual oxidase is presumably down-regulated to avoid host tissue damage by free reactive oxygen species during tissue restructuring. In contrast, the high pirk mRNA levels of infected compared to controls led us suspect an Imd signalling pathway inhibition [42]. Toll and Imd pathways were shown to be selectively activated depending on the composition of the pathogen-associated molecular patterns (PAMPs) with Toll pathway directed against Gram-positive bacteria and fungi and Imd pathway against Gram-negative bacterial cells [5]. Although, this mechanism of strict specificity was challenged, it may explain a reduced activity of the Imd pathway in a late response whereas Toll pathway, potentially selectively recognising the Gram-positive M. plutonius, seems to be sustained. Furthermore, activation of these pathways may differ in timing or be sequential [43, 44]. In other insects like Drosophila melanogaster, the Toll and Imd pathways are both known to be involved in the activation of cellular defence and melanisation response or at least overlap in regulation [44]. Therefore, this may also hold true for the honey bee, where we find a regulation of genes associated exactly with these defences.

Immune effector expression

Contrary to our expectations, we did not detect a significant differential immune response, specifically AMPs, in young infected larvae (3 days post infection = 72 h). In previous experiments, we showed that in 1 day-old and 3 days-old larvae infected per os with M. plutonius, an expressional change for a set of canonical immune genes of signal or effector proteins was not detectable [16]. Here, our transcriptomic analysis confirms these previous findings, although in future examinations of transcriptional AMP profiles a higher resolution of expressional timing would be desirable, particularly close to the infection time point. In comparison, P. larvae infections, causing AFB, the most dreaded bacterial disease of honey bee brood, elicit different immune effectors like AMPs 72 h after inoculation with the infecting agent [5, 40]. Although, we did not detect a significant change of ppo expression in infected larvae two putatively PPO cleaving serine proteases were up-regulated in 3 dpi M. plutonius infected larvae, presumably increasing the potential for active phenoloxidase (PO) in infection. A proteomic study on AFB diseased larvae demonstrated that in general, PO activity is highest on day 4 post-hatching (corresponding to 3 dpi in our study) and PPO seems to be accumulated in the haemolymph [21]. Although, PO is a well-established immune effector there may be also other functions as well as different pathways of activation. Moreover, several complement component associated protein genes were significantly altered in mRNA abundance with infection, including several C-type lectins (CTL1, CTL4, CTL6, CTL12) and another complement-like protein (LOC725964) in 3 dpi larvae, which have a potential role in opsonization of pathogens and possibly contributing to a PPO activation, a melanization/nodulation reaction and haemocyte-associated cellular defence [45].

On 6 dpi several immune effector genes were significantly up-regulated in infected individuals compared to controls, particularly hymenoptaecin and abaecin as well as two lysozymes. We did not detect defensin expressional regulation (Fig. 4B). Although, defensins usually are strongly up-regulated in adult bees upon infections [46] it seems to be less responsive in larvae. One reason may be the defensin content of the larval food, which is produced and secreted by the feeding workers [47]. For Drosophila melanogaster it was shown that both Toll and Imd pathways are necessary to up-regulate defensin [48]. In our set up, the Imd pathway seems to be inhibited and therefore prohibiting defensin expression is a mechanistic explanation for our observation. Moreover, defensins may also be expressed locally [49]. Hence, mRNA level variation in different tissues may remain undetected if measured in whole larvae like in our study.

By comparing the effect of time and treatment on immune effector mRNA levels, the time factor, and hence larval maturation, has by far a stronger effect on expression in both controls and infected individuals. This expressional regulation of immune effectors and other immune pathway genes is most likely due to the preparation for metamorphosis of the late-stage larvae like observed for other holometabolous insects [50]. This rational can also lead to an alternative interpretation of the observed gene expression differences. Hence, difference in immune effector abundance, particularly of AMPs, may also be attributed to a developmental retardation and a shift of developmental phases rather than a strong induced response to infection. This is supported by the expressional profiles of other gene groups like mentioned above. Unfortunately, we are not able to assess the strength of individual factors in this study.

However, the most striking transcriptional change in older infected larvae (both just M. plutonius and secondarily infected) was an increase in expression of cellular defence associated genes eater-like (LOC724421) and a scavenger receptor class B member 1 (LOC413408; also peste homologue) compared to controls. Eater is a haemocytes-specific Nimrod family cell surface receptor, in particular expressed by the plasmatocyte lineage in Drosophila melanogaster [51], a cell type analogous to honey bee granulocytes. Richardson et al. [52] reported that granulocytes are the dominant haemocyte type in honey bee larvae, so elevated expression of eater points to an increase in phagocytosis-associated cellular defence against Gram-positive bacteria. Moreover, Eater is the most crucial receptor in phagocytosis of E. faecalis, an Enterococcus species closely related to M. plutonius. The scavenger receptor class B member 1 is another protein gene that was shown to be involved in cell surface recognition and phagocytosis of pathogens [53].

Melanization cascade, complement-like and haemocyte-specific receptor gene expression as well as inhibition of the Imd pathway signalling genes integrate well in a mechanism with an activation of a cellular defence involved response in the larval haemolymph and possibly the gut. In Drosophila it was proposed that gut tissue residing haemocytes, which are spanning the gut epithelium, may be able to signal a gut infestation to the haemocoel and elicit a systemic response [48]. The honey bee larva is usually not able to stop multiplication of M. plutonius in the gut, but it still may be well prepared for a dissemination of bacterial cells into the haemolymph. Although, M. plutonius infections are most likely restricted to the larval gut and a spread of M. plutonius cells into the haemocoel was not shown so far, this phenomenon was already studied for other bacterial diseases like AFB [23]. An active breaching of the gut epithelium like demonstrated for P. larvae is not reported for M. plutonius, but a possible leakage of bacterial cells through loss of barrier function in late-stage infections could also be a plausible mechanism for EFB disease progression. Previously, loss of intestine barrier function was connected to enhanced immune system activation, inflammation, aging and reduced survival in Drosophila [54, 55]. While younger honey bee larvae seem to rely rather on PPO activation and opsonization of intruders like stated above, older larvae usually have higher numbers of haemocytes and therefore most likely invest in a stronger cellular defence [52]. This may contribute to an explanation why larvae survive better when infected in a later stage of development besides the reduced available time for multiplication of the pathogen inside the gut.

Behavioural change and hormone associated genes

Infections can have a significant effect on the behaviour of the host. Systemic changes during infection like shifts in metabolism and feeding behaviour are regulated by a network of hormones and neuropeptides, which usually also play a fundamental role in developmental control [56]. We found several indicators of expressional regulation of genes, which correlate with starvation and behavioural change. Several peptides, which are produced in specialized cells in the brain and peripheral tissues, are known to regulate hormone release and elicit behavioural change in insects [56]. Neuropeptide F was higher expressed in older infected larvae in both just primarily and secondarily infected larvae compared to controls while the serine protease gene scarface was down-regulated in these groups. This expressional pattern is usually associated with modification of feeding including a prolonged larval food uptake before metamorphosis, modulation of larval food foraging behaviour and may delay pupation like shown for other insects (e.g. Drosophila melanogaster) [30, 31]. Development and survival of larvae are also affected by the alteration of the synthesis of two hormones, JH and 20-HE, which play a crucial role in larval development and metamorphosis. Hence, differential expression of farnesol dehydrogenases genes as well as the down-regulated cytochrome P450 314A1 (also called shade) gene expression may point to an inhibition of the hydroxylation of ecdysone to active 20-HE and an altered JH synthesis. These steps of synthesis and modification are critical prior to the secretion and activity of these hormones. In our study, this is another indication of delay and an explanation for the regularly observed failure of developmental in EFB diseased larvae. An altered juvenile hormone as well as ecdysone synthesis gene expression are in concordance with this interpretation and underpin this view.

Secondary infection

Besides the primary infection with M. plutonius, in EFB several secondary invaders, i.e. opportunistic bacteria joining the causative agent in a later stage of infection, have been identified [2]. Previous studies obtained ambiguous results concerning an effect of the secondary invader, Paenibacillus alvei, on larval survival. Some approaches did not find increased mortality in in vitro infection experiments [7, 8] while Giersch et al. [57] could show an effect. Aside from an impact of secondary invaders on host survival there is a potential for an expressional response in secondary infection that does not directly translate to an apparently observable phenotype. In the present study, we did not find significantly increased mortality but a marked difference of expressed genes in secondary infection compared to only primary infected larvae, although a large proportion of DEGs overlapped. Interestingly, while most candidate genes are regulated in the same fashion in both contrasts, the secondary infection induces no significant effect or even a down-regulation of AMPs, some lysozymes and PRRs. Further, in contrast to solely M. plutonuis infected larvae lipoprotein genes like vg and A4, which are expected to have a function in immunity, were up-regulated. A secondary infection with P. alvei is apparently triggering a much stronger effect on metabolism by transcriptionally up-regulating a broad spectrum of different metabolic pathways compared to solely M. plutonius infected individuals. This response is potentially a specific tuning of the immune system and metabolism to a second invading bacterium in infected larvae.

Most noteworthy, significant enrichment of down-regulated hippo signalling pathway genes is indicating an interference of secondary infection with growth and development [58]. Unfortunately, in reference to several not or down-regulated canonical immune genes we cannot differentiate between a specific adjusted larval immune response to a secondary infection with a specific bacterium or a possible immune system inhibition by the secondary invader. This question has to be examined in additional experimental work.


The results of our analysis point to a change of the metabolic state and a dysregulation or temporal shift of structural and functional protein expression in EFB-associated bacterial infection of honey bee larvae. This transcriptional pattern supports the idea of an impairment of larval development as one of the key hallmarks of EFB disease. Although, we are aware of the potential involvement of toxins in EFB disease produced by several strains, we are confident that the illustration of infection we present here is the general underlying mechanism irrespective of toxin release by the pathogen. Late-stage infected larvae showed an increase of few canonical immune genes specifically antimicrobial peptides/lysozymes and haemocytic phagocytosis-associated genes. In contrast, in an earlier stage of infection no differential transcriptional immune effector response could be detected but rather an up-regulation of complement-like and melanization cascade gene expression. However, our results do not exclude potential immune reactivity at early stages of infection. Whether these changes and their timing are due to an adaptive defence mechanism, a manipulation or damage by the pathogen, or a side effect of either both, should be an aim for future experimental work. Several putatively regulatory non-coding RNAs may act as a part of the immune effector repertoire, signaling in infected individuals or be a result of physiological change and regulation but this requires further investigations. Further, we showed an up-regulation of specific immune effectors and cellular defence markers in preparation to metamorphosis, which was not shown before in the honey bee. An interference of infection with consequential immune system activation was rarely addressed in previous studies on transcriptional response of infected honey bee larvae. Therefore, our study provides new insights into this complex relation but more experimental examinations are needed to elucidate this interrelation and to isolate the involved factors.


Bacteria cultivation and larvae infection

M. plutonius strain 49.3 (wild type strain isolated in Switzerland, [8, 10] was cultivated and used for infection as described previously (OD600 adjusted to 0.3 with 2.4 × 107 CFU/ml; approx. 6.5 × 103 CFU per larva) [7]. Freshly hatched larvae were infected with M. plutonius and controls received sterile medium. On day 3 post infection a subset of larvae received a diet mixed with a secondary invader (P. alvei, [1]) or sterile medium (19:1). Paenibacillus alvei (LMG 13253) was provided by BCCM/LMG Bacteria Collection (Ghent University, Ghent, Belgium) and grown in specific medium as described in [59]. In a 12 day course of in vitro rearing, infected (only M. plutonius and M. plutonius + P. alvei) larvae had a significantly lower survival and larval weight compared to control groups (for more details see [7]).

Nucleic acid extraction and sample assessment

In total 3 to 5 replicate runs (24 larvae/replicate) per treatment were conducted. Samples were collected on day 3 post infection (3 dpi; n = 6 controls, n = 6 infected) and day 6 post infection (6 dpi; n = 6 controls, n = 6 infected, n = 6 secondary infected), respectively (Table S 1), weighed, sacrificed in liquid nitrogen and stored at − 80 °C until further processing. Total RNA and DNA were successively isolated from individual larvae. Sampled individuals were homogenized in 300 μl (3 dpi) and 600 μl (6 dpi) extraction buffer (NucleoSpin RNA purification kit, Macherey-Nagel, Düren, Germany), respectively. 50 μl (3 dpi) and 100 μl (6 dpi) of the homogenate was used for DNA extraction with the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). The remaining homogenate was used for RNA isolation using a NucleoSpin RNA purification kit (Macherey-Nagel, Düren, Germany) according to manufacturer’s protocol.

The pathogen load was estimated from isolated DNA samples with M. plutonius specific primers (sod A gene, forward: 5′-CAGCTAGTCGGTTTGGTTCC-3′; reverse: 5′-TTGGCTGTAGATAGAATTGACAAT-3′) relative to a honey bee reference gene (amCOI, forward: 5′-CCCCAGGATCATGAATTAGCAATGA-3′; reverse: 5′-TTCGGGGGAATGCTATATCAGGT-3′) via qPCR like described in [7], using M. plutonius genomic DNA as positive control as well as respective negative controls. Further, we determined pathogen loads as plasmid copy equivalents relative to a plasmid standard (provided by Eva Forsgren, SLU, Sweden) and took larval weight into account. Finally, we selected 4 larvae per treatment and time point for RNA sequencing, preferably with the highest RNA purity, determined on a NanoDrop 1000, and a similar pathogen load in the infected groups (or no detectable M. plutonius load for controls) (for further details on individual samples see Table S 9, Additional file 11).


The integrity of isolated RNA was verified (with the respective RIN) on an Agilent Bioanalyzer 2100 using an Agilent RNA 6000 Nano Kit as recommended by the manufacturer (Agilent Technologies, Waldbronn, Germany). The remaining genomic DNA was digested with TURBO DNase (Invitrogen, ThermoFischer Scientific, Paisley, United Kingdom). The amount of bacterial-derived rRNA sequences was reduced with the Ribo-Zero magnetic kit (Epicentre Biotechnologies, Madison, WI, USA) and the removal of eukaryotic-derived rRNA sequences was performed with the NEBNext Poly(A) mRNA Magnetic Isolation Module (New England BioLabs, Frankfurt am Main, Germany). For sequencing, cDNA libraries were constructed with a NEBNext Ultra II RNA Library Prep Kit for Illumina (New England BioLabs, Frankfurt am Main, Germany). Library sizes and quality of samples were assessed on an Agilent Bioanalyzer 2100 using an Agilent High Sensitivity DNA Kit as recommended by the manufacturer (Agilent Technologies, Waldbronn, Germany). Determination of the concentration of libraries was performed on the Qubit® dsDNA HS Assay Kit as recommended by the manufacturer (Life Technologies GmbH, Darmstadt, Germany). Sequencing was performed on a HiSeq2500 instrument (Illumina Inc., San Diego, CA, USA) using the HiSeq Rapid PE Cluster Kit v2 for cluster generation and the HiSeq Raid SBS Kit (300 cycles) for sequencing in the paired-end mode (read length of 150 bp) and running 2 × 150 cycles.

Bacterial community analysis

Bacterial 16S rRNA gene amplicons were generated using fusion primers TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG (MiSeq_B_V3_for_klindworth) and GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC (MiSeq_B_V4_rev_klindworth) including bacteria targeting primers from Klindworth et al. [60]. The PCR reaction mixture with a total volume 50 μl contained 1 U Phusion high fidelity DNA polymerase (Biozym Scientific, Oldendorf, Germany), 5% DMSO, 0.2 mM of each primer, 200 μM dNTPs, 0.2 μl of 50 mM MgCl2, and 25 ng of isolated DNA. Thermal cycling scheme for bacterial amplicons was as follows: initial denaturation for 1 min at 98 °C, 25 cycles at 98 °C for 45 s, 45 s at 60 °C, and 30 s at 72 °C, and a final extension at 72 °C for 5 min. The resulting PCR products were checked by agarose gel electrophoresis for appropriate size and purified using the MagSi-NGSPREP Plus Magnetic beads (Steinbrenner Laborsysteme GmbH, Wiesenbach, Germany) as recommended by the manufacturer. PCR products were used to attach indices and Illumina sequencing adapters using the Nextera XT Index kit (Illumina, San Diego). Index PCR was performed using 5 μl of template PCR product, 2.5 μl of each index primer, 12.5 μl of 2× KAPA HiFi HotStart ReadyMix and 2.5 μl PCR grade water. Thermal cycling scheme was as follows: 95 °C for 3 min, 8 cycles of 30 s at 95 °C, 30 s at 55 °C and 30 s at 72 °C and a final extension at 72 °C for 5 min. Quantification of the products was performed using the Quant-iT dsDNA HS assay kit and a Qubit fluorometer (Invitrogen GmbH, Karlsruhe, Germany) following the manufacturer’s instructions. MagSi-NGSPREP Plus Magnetic beads (Steinbrenner Laborsysteme GmbH, Wiesenbach, Germany) were used for purification of the indexed products as recommended by the manufacturer and normalization was performed using the Janus Automated Workstation from Perkin Elmer (Perkin Elmer, Waltham Massachusetts, USA). Sequencing was conducted using the Illumina MiSeq platform with dual indexing and MiSeq reagent kit v3 (600 cycles) as recommended by the manufacturer (Illumina). CASAVA software (Illumina) was used for demultiplexing and clipping of adapter sequences from the raw amplicon sequences. Fastp (ver. 0.19.6) [61] with a minimum phred score of 20, a minimum length of 50 base pairs, a sliding window size of four bases, was used for read correction by overlap and adapter removal of the Illumina Nextera primers. Paired-end reverse reads were merged with the paired end read merger (PEAR ver. 0.9.11) [62] with default settings. Additionally, reverse and forward primer sequences were removed with cutadapt (ver. 1.18) [63] with default settings. Sequences were then size filtered (≤300 bp were removed) and dereplicated by vsearch (ver. 2.11.) [64]. Denoising was performed with the UNOISE3 module of vsearch and a set minimum size of 8 reads. Chimeric sequences were excluded with the UCHIME module of vsearch. This included de novo chimera and reference-based chimera removal against the SILVA SSU 132 NR database [65, 66]. Taxonomy assignments were performed with BLASTn (ver. 2.9.0+) [67] against the SILVA SSU 138 NR database with an identity threshold of 90%. Abundance bar charts were created from amplicon sequence variants with the ggplot2 package with R (R core team, 2017) and RStudio® [68].

The operational taxonomic unit (OTU) analysis confirmed that more than 99% of bacteria in the infected individuals were represented by M. plutonius reads compared to the control larvae (free of M. plutonius) in both 3 days and 6 days post infection individuals. We determined the relative compositions of the larval microbiota, which was composed of typical honey bee specific (e.g. Fructobacillus, Lactobacillus, Bombella, etc.) and also environmental bacterial genera (Fig. S 5, Additional file 1).

Transcriptome analysis

We performed our analysis with trimmed paired-end raw reads. Mapping and transcript abundance quantification was carried out with Salmon (ver. 2.1.0) [69] employing the most recent available Apis mellifera reference genome (Amel_HAv3.1, Annotation release 104). We recovered a high read number (total fragments: 14.89 × 106 ± 9.0%) and high mapping efficiency (mapping rate range: 77.3–84.6%), which was consistent overall samples (Table S 9, Additional file 11) and complies with common quality requirements. The analysis of the quantified transcripts was followed up with the ThreeDRNAseq R package (ver. 2.0.1) [70]. Read counts and transcript per million reads (TPM) were generated using tximport R package version 1.10.0 and lengthScaledTPM method [71] with inputs of transcript quantifications from Salmon. Overall, we found 27,774 expressed transcripts of 13,829 genes. Low expressed transcripts and genes were filtered based on analysis of the data mean-variance trend (count per million reads (CPM) cut-off ≥1 and min. Sample cut-off ≥4; 4 samples resembling one group, so the chance that a very low expressed transcript is most likely to be significant if expressed in the same group with low variation) yielding 17,751 transcripts of 10,394 genes (Fig. S 6 and Fig. S 7, Additional file 1). Gene and transcript read counts were normalised to log2-CPM using the TMM method (Fig. S 8, Additional file 1; for complete read count and TPM data see Table S 10, Additional file 13) [72]. Further, in-group consistency and between-group variation was assessed with a principal component analysis (Fig. 1) and batch-effects were corrected with RUVr function from the RUV R package (ver. [73]. The contrast groups in Table 1 were used to compare the expression changes between conditions of experimental design. For DEGs, the log2 fold change of transcript abundance were calculated based on contrast groups and significance of expression changes were determined using t-test. P-values of multiple testing were adjusted with Benjamini-Hochberg method to correct false discovery rate (FDR). A DEG was significant in a contrast group if it had adjusted p-value < 0.05 and log2 FC > 1 (for full account see Table S 11, Additional file 13).

The functional annotation clustering analysis was performed with DAVID 6.8 [74]. DEGs were converted to the respective DAVID data base accessions using the gene accession conversion tool. GO term grouping (molecular function, biological process, cellular compartment), the KEGG biological pathways as well as UP Keywords and InterPro categories were included in clustering and the threshold of significant terms was set to p < 0.1 (for further settings see Table S 2, Additional file 4). The results of the FAC were processed and visualized with the EnrichmentMap (ver. 2.0) plugin for Cytoscape (ver. 3.8) to illustrate and contrast the cluster and term associations of the compared groups. Functional terms analysis of secondary infection (Cont vs. Inf and Cont vs. secInf) was carried out with sets of input DEGs under relaxed conditions (FDR < 0.05; − 0.585 < log2FC > 0.585; Table S 12, Additional file 14), up- and down-regulated genes separated, to increase resolution and then proceeded like described above.

Correlation analysis

In M. plutonius infected larvae, weight is an important factor and indicator of developmental progression and infection status. We realized that although our M. plutonius specific primers were reliable in detecting bacterial DNA, we were not able to differentiate whether bacterial nucleic acids originated from live bacteria or accumulated from dead cells. Weight and bacterial load did not correlate well in our full set of infected individuals (Pearson r = 0.15, p = 0.62), although the subset of solely M. plutonius infected correlated stronger (Pearson r = 0.66, p = 0.08; Fig. S 1, Additional file 1). Hence, we did not include infection status and performed an independent correlation analysis of mRNA levels and larval weights to avoid misinterpretation of infection loads and weight changes. Moreover, this procedure enabled us to analyse infection load independent expressional changes in infected larvae. Significance thresholds of Pearson correlations of expressed gene TPMs and weight were adjusted for multiple testing (Benjamin-Hochberg: FDR < 0.05) and the functional annotation tool in the DAVID database was used to retrieve information on functions of expressed genes significantly correlated with weight with the same procedure as before (see above).

Availability of data and materials

The dataset supporting the findings of this article is available in the NCBI SRA database under the BioProject accession number PRJNA745504 (





American foulbrood


Antimicrobial peptide


Chitin-binding protein


Counts per million


Cellular retinaldehyde-binding protein and triple functional protein domain


Differentially expressed genes


Days post infection


Domain of unknown function 1676


extracellular matrix


European foulbrood


Functional annotation clustering


Fold change


False discovery rate


Gene ontology


Juvenile hormone


Kyoto encyclopedia of genes and genomes


Operational taxonomic unit


Pathogen-associated molecular pattern


Principal component analysis


Peritrophic matrix


Phenoloxidase/ prophenoloxidase


Pattern recognition receptor


RNA integrity number


Transcripts per million




Very high density lipoprotein


Zona pellucida


  1. Forsgren E. European foulbrood in honey bees. J Invertebr Pathol. 2010;103(Suppl 1):S5–9.

    Article  Google Scholar 

  2. Forsgren E, Locke B, Sircoulomb F, Schäfer MO. Bacterial diseases in honeybees. Curr Clin Micro Rpt. 2018;5:18–25.

    Article  Google Scholar 

  3. Potts Z, Imperatriz-Fonseca V, Ngo H, Aizen MA, Biesmeijer JC, Breeze TD, et al. Safeguarding pollinators and their values to human well-being. Nature. 2016;540:220–9.

    Article  CAS  Google Scholar 

  4. Mutinelli F. The spread of pathogens through trade in honey bees and their products (including queen bees and semen): overview and recent developments. Rev Sci Tech. 2011;30:257–71.

    Article  CAS  Google Scholar 

  5. Evans JD, Aronstein K, Chen YP, Hetru C, Imler J-L, Jiang H, et al. Immune pathways and defence mechanisms in honey bees Apis mellifera. Insect Mol Biol. 2006;15:645–56.

    Article  CAS  Google Scholar 

  6. Kathe E, Seidelmann K, Lewkowski O, Le Conte Y, Erler S. Changes in chemical cues of Melissococcus plutonius infected honey bee larvae. Chemoecology. 2021;31:189–200.

    Article  CAS  Google Scholar 

  7. Lewkowski O, Erler S. Virulence of Melissococcus plutonius and secondary invaders associated with European foulbrood disease of the honey bee. MicrobiologyOpen. 2019;8:e00649.

    Article  CAS  Google Scholar 

  8. Grossar D, Kilchenmann V, Forsgren E, Charrière J-D, Gauthier L, Chapuisat M, et al. Putative determinants of virulence in Melissococcus plutonius, the bacterial agent causing European foulbrood in honey bees. Virulence. 2020;11:554–67.

    Article  CAS  Google Scholar 

  9. Budge GE, Barrett B, Jones B, Pietravalle S, Marris G, Chantawannakul P, et al. The occurrence of Melissococcus plutonius in healthy colonies of Apis mellifera and the efficacy of European foulbrood control measures. J Invertebr Pathol. 2010;105:164–70.

    Article  Google Scholar 

  10. Djukic M, Erler S, Leimbach A, Grossar D, Charrière J-D, Gauthier L, et al. Comparative genomics and description of putative virulence factors of Melissococcus plutonius, the causative agent of European foulbrood disease in honey bees. Genes. 2018;9:419.

    Article  CAS  Google Scholar 

  11. Sturtevant AP. A study of the behavior of bees in colonies affected by European foulbrood. Washington, D.C: U.S. Deptartment of Agriculture; 1920.

    Book  Google Scholar 

  12. White GF. European foulbrood. U.S. Department of Agriculture, bulletin no. 810. Washington: Government Printing Office; 1920.

    Google Scholar 

  13. Tarr HLA. Studies on European foul brood of bees II. The production of the disease experimentally. Ann Appl Biol. 1936;23:558–84.

    Article  Google Scholar 

  14. Bailey L. European foulbrood. Am Bee J. 1961;101:89–92.

    Google Scholar 

  15. Schmid-Hempel P. Evolutionary parasitology: the integrated study of infections, immunology, ecology, and genetics. Oxford: Oxford University Press; 2011.

    Book  Google Scholar 

  16. Ehrenberg S, Lewkowski O, Erler S. Dyeing but not dying: Colourful dyes as a non-lethal method of food labelling for in vitro-reared honey bee (Apis mellifera) larvae. J Insect Physiol. 2019;113:1–8.

    Article  CAS  Google Scholar 

  17. Takamatsu D, Osawa A, Nakamura K, Yoshiyama M, Okura M. High-level resistance of Melissococcus plutonius clonal complex 3 strains to antimicrobial activity of royal jelly. Environ Microbiol Rep. 2017;9:562–70.

    Article  CAS  Google Scholar 

  18. Bailey L. Melissococcus pluton, the cause of European foulbrood of honey bees (Apis spp.). J Appl Bacteriol. 1983;55:65–9.

    Article  Google Scholar 

  19. Tarr HLA. Studies on European foul brood of bees IV. On the attempted cultivation of Bacillus pluton, the susceptibility of individual larvae to inoculation with this organism and its localization within its host. Ann Appl Biol. 1938;25:815–21.

    Article  Google Scholar 

  20. Takamatsu D, Sato M, Yoshiyama M. Infection of Melissococcus plutonius clonal complex 12 strain in European honeybee larvae is essentially confined to the digestive tract. J Vet Med Sci. 2016;78:29–34.

    Article  CAS  Google Scholar 

  21. Chan QW, Melathopoulos AP, Pernal SF, Foster LJ. The innate immune and systemic response in honey bees to a bacterial pathogen, Paenibacillus larvae. BMC Genomics. 2009;10:387.

    Article  CAS  Google Scholar 

  22. Nakamura K, Okumura K, Harada M, Okamoto M, Okura M, Takamatsu D. Peritrophic matrix-degrading proteins are dispensable virulence factors in a virulent Melissococcus plutonius strain. Sci Rep. 2021;11:8798.

    Article  CAS  Google Scholar 

  23. Djukic M, Brzuszkiewicz E, Fünfhaus A, Voss J, Gollnow K, Poppinga L, et al. How to kill the honey bee larva: genomic potential and virulence mechanisms of Paenibacillus larvae. PLoS One. 2014;9:e90914.

    Article  CAS  Google Scholar 

  24. Danty E, Arnold G, Burmester T, Huet JC, Huet D, Pernollet JC, et al. Identification and developmental profiles of hexamerins in antenna and hemolymph of the honeybee, Apis mellifera. Insect Biochem Molec. 1998;28:387–97.

    Article  CAS  Google Scholar 

  25. Shipman BA, Ryan RO, Schmidt JO, Law JH. Purification and properties of a very high density lipoprotein from the hemolymph of the honeybee Apis mellifera. Biochemistry. 1987;26:1885–9.

    Article  CAS  Google Scholar 

  26. Sutherland TD, Campbell PM, Weisman S, Trueman HE, Sriskantha A, Wanjura WJ, et al. A highly divergent gene cluster in honey bees encodes a novel silk family. Genome Res. 2006;16:1414–21.

    Article  CAS  Google Scholar 

  27. Smith CR, Morandin C, Noureddine M, Pant S. Conserved roles of Osiris genes in insect development, polymorphism and protection. J Evol Biol. 2018;31:516–29.

    Article  CAS  Google Scholar 

  28. Jovine L, Darie CC, Litscher ES, Wassarman PM. Zona pellucida domain proteins. Annu Rev Biochem. 2005;74:83–114.

    Article  CAS  Google Scholar 

  29. Lee SY, Kwon TH, Hyun JH, Choi JS, Kawabata SI, Iwanaga S, et al. In vitro activation of pro-phenol-oxidase by two kinds of pro-phenol-oxidase-activating factors isolated from hemolymph of coleopteran, Holotrichia diomphalia larvae. Eur J Biochem. 1998;254:50–7.

    Article  CAS  Google Scholar 

  30. Wu Q, Wen T, Lee G, Park JH, Cai HN, Shen P. Developmental control of foraging and social behavior by the Drosophila neuropeptide Y-like system. Neuron. 2003;39:147–61.

    Article  CAS  Google Scholar 

  31. Prasad N, Hens K. Sugar promotes feeding in flies via the serine protease homolog scarface. Cell Rep. 2018;24:3194–206.e4.

    Article  CAS  Google Scholar 

  32. Bomtorin AD, Mackert A, Rosa GC, Moda LM, Martins JR, Bitondi MM, et al. Juvenile hormone biosynthesis gene expression in the corpora allata of honey bee (Apis mellifera L.) female castes. PLoS One. 2014;9(1):e86923.

    Article  CAS  Google Scholar 

  33. Morris K, Mattick J. The rise of regulatory RNA. Nat Rev Genet. 2014;15:423–37.

    Article  CAS  Google Scholar 

  34. Zhang Y, Cao X. Long noncoding RNAs in innate immunity. Cell Mol Immunol. 2016;13:138–47.

    Article  Google Scholar 

  35. Chan QW, Foster LJ. Changes in protein expression during honey bee larval development. Genome Biol. 2008;9:R156.

    Article  CAS  Google Scholar 

  36. Merzendorfer H, Zimoch L. Chitin metabolism in insects: structure, function and regulation of chitin synthases and chitinases. J Exp Biol. 2003;206:4393–412.

    Article  CAS  Google Scholar 

  37. Hynes RO. The extracellular matrix: not just pretty fibrils. Science. 2009;326:1216–9.

    Article  CAS  Google Scholar 

  38. Santos DE, Alberici LC, Hartfelder K. Mitochondrial structure and dynamics as critical factors in honey bee (Apis mellifera L.) caste development. Insect Biochem Mol Biol. 2016;73:1–11.

    Article  CAS  Google Scholar 

  39. Arrese EL, Soulages JL. Insect fat body: energy, metabolism, and regulation. Annu Rev Entomol. 2010;55:207–25.

    Article  CAS  Google Scholar 

  40. Cornman RS, Lopez D, Evans JD. Transcriptional response of honey bee larvae infected with the bacterial pathogen Paenibacillus larvae. PLoS One. 2013;8:e65424.

    Article  CAS  Google Scholar 

  41. Haslbeck M, Weinkauf S, Buchner J. Small heat shock proteins: simplicity meets complexity. J Biol Chem. 2019;294:2121–32.

    Article  CAS  Google Scholar 

  42. Kleino A, Silverman N. Regulation of the Drosophila Imd pathway by signaling amyloids. Insect Biochem Mol Biol. 2019;108:16–23.

    Article  CAS  Google Scholar 

  43. Tanji T, Hu X, Weber AN, Ip YT. Toll and IMD pathways synergistically activate an innate immune response in Drosophila melanogaster. Mol Cell Biol. 2007;27:4578–88.

    Article  CAS  Google Scholar 

  44. De Gregorio E, Spellman PT, Tzou P, Rubin GM, Lemaitre B. The toll and Imd pathways are the major regulators of the immune response in Drosophila. EMBO J. 2002;21:2568–79.

    Article  Google Scholar 

  45. Xia X, You M, Rao XJ, Yu XQ. Insect C-type lectins in innate immunity. Dev Comp Immunol. 2018;83:70–9.

    Article  CAS  Google Scholar 

  46. Ilyasov RA, Gaifullina LR, Saltykova ES, Poskryakov AV, Nikolaenko AG. Defensins in the honeybee antiinfectious protection. J Evol Biochem Phys. 2013;49:1–9.

    Article  CAS  Google Scholar 

  47. Fujiwara S, Imai J, Fujiwara M, Yaeshima T, Kawashima T, Kobayashi K. A potent antibacterial protein in royal jelly. Purification and determination of the primary structure of royalisin. J Biol Chem. 1990;265:11333–7.

    Article  CAS  Google Scholar 

  48. Charroux B, Royet J. Elimination of plasmatocytes by targeted apoptosis reveals their role in multiple aspects of the Drosophila immune response. Proc Natl Acad Sci U S A. 2009;106:9797–802.

    Article  Google Scholar 

  49. Tzou P, Ohresser S, Ferrandon D, Capovilla M, Reichhart JM, Lemaitre B, et al. Tissue-specific inducible expression of antimicrobial peptide genes in Drosophila surface epithelia. Immunity. 2000;13:737–48.

    Article  CAS  Google Scholar 

  50. Johnston PR, Paris V, Rolff J. Immune gene regulation in the gut during metamorphosis in a holo- versus a hemimetabolous insect. Phil Trans R Soc B. 2019;374:20190073.

    Article  CAS  Google Scholar 

  51. Kocks C, Cho JH, Nehme N, Ulvila J, Pearson AM, Meister M, et al. Eater, a transmembrane protein mediating phagocytosis of bacterial pathogens in Drosophila. Cell. 2005;123:335–46.

    Article  CAS  Google Scholar 

  52. Richardson RT, Ballinger MN, Qian F, Christman JW, Johnson RM. Morphological and functional characterization of honey bee, Apis mellifera, hemocyte cell communities. Apidologie. 2018;49:397–410.

    Article  CAS  Google Scholar 

  53. Philips JA, Rubin EJ, Perrimon N. Drosophila RNAi screen reveals CD36 family member required for mycobacterial infection. Science. 2005;309:1251–3.

    Article  CAS  Google Scholar 

  54. Rera M, Clark RI, Walker DW. Intestinal barrier dysfunction links metabolic and inflammatory markers of aging to death in Drosophila. Proc Natl Acad Sci U S A. 2012;109:21528–33.

    Article  Google Scholar 

  55. Rera M, Azizi MJ, Walker DW. Organ-specific mediation of lifespan extension: more than a gut feeling? Ageing Res Rev. 2013;12:436–44.

    Article  CAS  Google Scholar 

  56. Pool AH, Scott K. Feeding regulation in Drosophila. Curr Opin Neurobiol. 2014;29:57–63.

    Article  CAS  Google Scholar 

  57. Giersch T, Barchia I, Hornitzky M. Can fatty acids and oxytetracycline protect artificially raised larvae from developing European foulbrood? Apidologie. 2010;41:151–9.

    Article  CAS  Google Scholar 

  58. Ashby R, Forêt S, Searle I, Maleszka R. MicroRNAs in honey bee caste determination. Sci Rep. 2016;6:18794.

    Article  CAS  Google Scholar 

  59. Erler S, Denner A, Bobiş O, Forsgren E, Moritz RFA. Diversity of honey stores and their impact on pathogenic bacteria of the honeybee, Apis mellifera. Ecol Evol. 2014;4:3960–7.

    Article  Google Scholar 

  60. Klindworth A, Pruesse E, Schweer T, Peplies J, Quast C, Horn M, et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 2013;41:e1.

    Article  CAS  Google Scholar 

  61. Chen S, Zhou Y, Chen Y, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90.

    Article  CAS  Google Scholar 

  62. Zhang J, Kobert K, Flouri T, Stamatakis A. PEAR: a fast and accurate Illumina paired-end reAd mergeR. Bioinformatics. 2014;30:614–20.

    Article  CAS  Google Scholar 

  63. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17(10–2).

  64. Rognes T, Flouri T, Nichols B, Quince C, Mahé F. VSEARCH: a versatile open source tool for metagenomics. PeerJ. 2016;4:e2584.

    Article  Google Scholar 

  65. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet C, Ghalith GA, et al. QIIME 2: reproducible, interactive, scalable, and extensible microbiome data science. Peer J Preprints. 2018;6:e27295v2.

    Article  Google Scholar 

  66. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:D590–6.

    Article  CAS  Google Scholar 

  67. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.

    Article  CAS  Google Scholar 

  68. Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer-Verlag; 2016.

    Book  Google Scholar 

  69. Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14:417–9.

    Article  CAS  Google Scholar 

  70. Guo W, Tzioutziou NA, Stephen G, Milne I, Calixto CPG, Waugh R, et al. 3D RNA-seq: a powerful and flexible tool for rapid and accurate differential expression and alternative splicing analysis of RNA-seq data for biologists. RNA Biol. 2021;18:1574–87.

    Article  CAS  Google Scholar 

  71. Soneson C, Matthes KL, Nowicka M, Law CW, Robinson MD. Isoform prefiltering improves performance of count-based methods for analysis of differential transcript usage. Genome Biol. 2016;17:12.

    Article  CAS  Google Scholar 

  72. Bullard JH, Purdom E, Hansen KD, et al. Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatics. 2010;11:94.

    Article  CAS  Google Scholar 

  73. Risso D, Ngai J, Speed TP, Dudoit S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat Biotechnol. 2014;32:896–902.

    Article  CAS  Google Scholar 

  74. Huang D, Sherman B, Lempicki R. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.

    Article  CAS  Google Scholar 

Download references


We are grateful to Sandra Ehrenberg for help with lab work, as well as Jean-Daniel Charrière (Swiss Bee Research Center, Agroscope) for providing M. plutonius strain 49.3 and Eva Forsgren (Swedish University of Agricultural Sciences, Uppsala) for providing the plasmid used for M. plutonius quantification. Further, we would like to thank Dominik Schneider for providing the amplicon pipeline.


Open Access funding enabled and organized by Projekt DEAL. Financial support was granted by the German Research Foundation - DFG (ER 786/1–1 to S.E.).

Author information

Authors and Affiliations



OL & SE conceived the scientific questions and experimental design. OL collected larvae samples and extracted the RNA. AP performed the RNA quality analysis and sequencing. OL performed the transcriptomic data analysis, clustering analysis and wrote the original draft of the manuscript. SE led project administration and supervision. All authors contributed significantly to the writing of the manuscript and approved it for publication.

Corresponding authors

Correspondence to Oleg Lewkowski or Silvio Erler.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1 

: Fig. S1 - S10. For details see Figure captions.

Additional file 2

: Fig. S11. Detailed network figure of FAC analysis clusters and term associations.

Additional file 3

: Table S1. Correlation analysis of weight and mRNA abundance.

Additional file 4

: Table S2. Detailed results of FAC analysis for temporal DEGs in control (Cont 3 dpi vs 6 dpi) and infected (Inf 3 dpi vs 6 dpi) group, respectively.

Additional file 5

: Table S3. Expressional comparison of chitin-binding/cuticula-associated protein genes.

Additional file 6

: Table S4. Expressional comparison of DUF1676 domian containing genes and the Osiris family.

Additional file 7

: Table S5. Expressional comparison of genes related to the canonical immune pathways Toll and Imd/JNK.

Additional file 8

: Table S6. Expressional comparison of candidate genes with a function in A) silk formation, storage and/or transport of energy reserves and B) related to immune response corresponding to Fig. 4.

Additional file 9

: Table S7. Expressional comparison of genes of enzymes involved in insect hormone biosynthesis.

Additional file 10

: Table S8. Expressional comparison of heme domain-containing protein genes.

Additional file 11

: Table S9. Summary of read mapping for all samples performed with mapping-based mode and a transcriptome index in Salmon.

Additional file 12

: Table S10. Read counts and TPM reads. Raw counts were generated from transcript quantifications from Salmon using tximport R package and TPMs with lengthScaledTPM method.

Additional file 13

: Table S11. List of significant DEGs by contrast group category with adjusted p-value < 0.05 and log2 FC > 1.

Additional file 14

: Table S12. Detailed results of FAC analysis for DEGs of 6 dpi control (Cont) versus primarily infected (Inf) and secondary infected (secInf) group, respectively.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lewkowski, O., Poehlein, A., Daniel, R. et al. In the battle of the disease: a transcriptomic analysis of European foulbrood-diseased larvae of the Western honey bee (Apis mellifera). BMC Genomics 23, 837 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: