Infection routes matter in population-specific responses of the red flour beetle to the entomopathogen Bacillus thuringiensis
© Behrens et al.; licensee BioMed Central Ltd. 2014
Received: 14 January 2014
Accepted: 4 June 2014
Published: 7 June 2014
Pathogens can infect their hosts through different routes. For studying the consequences for host resistance, we here used the entomopathogen Bacillus thuringiensis and the red flour beetle Tribolium castaneum for oral and systemic (i. e. pricking the cuticle) experimental infection. In order to characterize the molecular mechanisms underpinning the two different infection routes, the transcriptomes of beetles of two different T. castaneum populations – one recently collected population (Cro1) and a commonly used laboratory strain (SB) – were analyzed using a next generation RNA sequencing approach.
The genetically more diverse population Cro1 showed a significantly larger number of differentially expressed genes. While both populations exhibited similar reactions to pricking, their expression patterns in response to oral infection differed remarkably. In particular, the Cro1 population showed a strong response of cuticular proteins and developmental genes, which might indicate an adaptive developmental flexibility that was lost in the SB population presumably as a result of inbreeding. The immune response of SB was primarily based on antimicrobial peptides, while Cro1 relied on responses mediated by phenoloxidase and reactive oxygen species, which may explain the higher resistance of this strain against oral infection.
Our data demonstrate that immunological and physiological processes underpinning the two different routes of infection are clearly distinct, and that host populations particularly differ in responses to oral infection. Furthermore, gene expression upon pricking infection entailed a strong signal of wounding, highlighting the importance of pricking controls in future infection studies.
KeywordsTribolium castaneum Bacillus thuringiensis Host-parasite interactions Oral infection Pricking infection RNA sequencing Transcriptome
The route by which pathogens infect their hosts can have important consequences for host-pathogen interactions . For example, a recent meta-analysis showed that virulence was higher in pathogens infecting wounded skin, compared with those ingested or inhaled . Pathogens that can infect through alternative routes are particularly interesting systems with which to test the evolutionary and physiological consequences of infection through different routes. A recent study investigated the evolution of resistance of Drosophila melanogaster hosts against Pseudomonas entomophila bacteria upon oral as compared to systemic (i. e. pricking the cuticle) infection . Interestingly, hosts evolved resistance towards the bacteria for both routes of infection. However, there was no cross-resistance, i. e. D. melanogaster selection lines that had evolved resistance against P. entomophila upon oral infection were not more resistant against P. entomophila upon pricking and vice versa. This route-specificity indicates that the physiological underpinnings of resistance and therefore the evolutionary trajectories of adaptation differ for the routes of infection. However, it is currently unclear whether such route-specificity is a general phenomenon or restricted to this particular host-pathogen system.
In the present study, we studied infection route-specificity in the red flour beetle Tribolium castaneum upon oral and pricking infection with the entomopathogen Bacillus thuringiensis[4, 5]. We compared gene expression between pricking and oral infection, using an Illumina next generation sequencing approach (RNA-seq). RNA-seq is a powerful tool that enables a very precise quantification of transcript levels on a genome-wide scale . For our comparative RNA-seq based study, a sterile wounding treatment was included to distinguish between effects of wounding alone and bacterial infection.
T. castaneum is a relevant global stored product pest that has developed into a fully-fledged insect model organism . The genome of the T. castaneum strain Georgia 2 (GA-2) has been fully sequenced in 2008 . B. thuringiensis is a Gram-positive bacterium that forms highly resistant endospores and plasmid-encoded crystalline inclusions (Cry proteins), which are toxic upon oral ingestion . Cry toxins provide specificity for insect orders upon oral infection . We recently showed that spore-crystal mixtures of the B. thuringiensis bv.tenebrionis strain are infectious to T. castaneum upon oral exposure . However, B. thuringiensis is also able to efficiently infect and kill T. castaneum upon pricking infection . Insects may regularly suffer from wounding in their natural environments, such that infection through the wounded cuticle is likely to occur in nature. Moreover, septic pricking in the laboratory can lead to spore production in cadavers (unpublished data). Defense of T. castaneum against septic infection with B. thuringiensis has been studied as a model for eco-immunology and the host has been shown to be capable of specific immune priming within and across generations [5, 11]. Recently, it was also shown that the oral exposure to B. thuringiensis spore supernatants leads to priming of T. castaneum larvae. The analysis of the host responses to oral and pricking infection may contribute to an enhanced understanding of the mechanistic underpinnings of this astonishing degree of immunological adaptation.
Natural variation in resistance may provide a basis for studying evolved differences among host populations, and thereby provides a potentially important source of information for the identification of the genetic causes of resistance. We have previously shown that a recently collected population of T. castaneum (Cro1) showed enhanced resistance to oral infection with B. thuringiensis, as compared to commonly used laboratory populations (SB, GA-2) . In the present study, we therefore compared gene expression between the populations Cro1 and SB upon oral and pricking infection with B. thuringiensis. We found a higher number of differentially expressed genes in the Cro1 population for both routes of infection. Intriguingly, gene expression profiles differed strongly between the oral and pricking routes of infection, indicating that immunological and physiological processes underpinning the two routes of infection are clearly distinct. Furthermore, we demonstrate that pricking without bacteria, i. e. a sterile aseptic wounding, leads to strong immune activation and therefore represents a necessary control for pricking infections.
Single-nucleotide resolution transcriptome of T. castaneumby RNA-seq
Numbers of differentially expressed genes differ between populations and infection routes
Cuticle genes are highly enriched among differentially expressed genes
Populations show specific patterns of gene regulation
In order to further investigate the enrichment of the GO term “structural constituent of cuticle” in sets of differentially expressed genes described in the preceding paragraph, we also applied the overrepresentation analysis to a list of cuticle genes, i. e. to a list of all genes with the GO term “structural constituent of cuticle”; see Additional file 8. For BttP, Figure 5 shows a similar pattern in the two different populations: cuticle genes are highly enriched 18 h after infection in Cro1 and SB in both significantly up- and downregulated genes. For PC, we observe an enrichment of cuticle genes only in the lists of significantly downregulated genes of the SB treatments. Noticeably, Figure 5 demonstrates that cuticle genes are highly enriched in the significantly upregulated genes 6 and 18 h after oral infection in Cro1 (and not in the downregulated genes) while for SB, they are significantly overrepresented in the significantly downregulated genes 6 and 18 h after oral infection (and not in the upregulated genes). This clearly reveals another population-specific reaction upon oral infection: for Cro1, relatively more cuticle genes are up- than downregulated and for SB, relatively more cuticle genes are down- than upregulated.
A recent transcriptome study of T. castaneum hypothesized a crosstalk between stress-related genes, cytochrome P450s and between immune genes . This is why we defined a set of “stress genes” (see Additional file 9) and another set of “P450 genes” (see Additional file 10) and tested for enrichment of these gene sets in our differential expression results. Figure 5 shows that stress genes are significantly overrepresented in the lists of significantly upregulated genes (5/12 lists) with a slightly stronger enrichment in Cro1 compared to SB. Thus, this supports the hypothesis of an interdependency between stress and immunity related genes. Remarkably, we observe a very strong overrepresentation of cytochrome P450s in the lists of significantly downregulated genes (9/12 lists). Hence, we suggest that downregulation of cytochrome P450 genes plays an important role in response to all three treatments BttO, BttP and PC. Interestingly, for Cro1, P450 genes are enriched in the 6 h as well as 18 h after exposure treatments, while for SB, P450 genes are only enriched in the 6 h after exposure treatments. Thus, we hypothesize a prolonged downregulation of cytochrome P450s in Cro1 compared to SB.
Upon oral infection, many odorant binding proteins are found differentially regulated (see Additional file 11). Odorant binding proteins have been associated to host defense and there is evidence that they play a role in resistance to B. thuringiensis intoxication in T. castaneum. For this reason a further set of “odorant binding genes” was defined; see Additional file 12. As can be observed in Figure 5, odorant binding genes are only enriched among significantly upregulated genes in SB, 18 h after pricking infection. Nevertheless, we found several odorant binding genes differentially expressed in BttO, BttP and PC. Moreover, each population expresses a different set of odorant binding genes. For example, 6 h after oral infection, OBP-C04, OBP-10 and OBP-C12 are significantly upregulated in Cro1 while OBP-8, OBP-12 and OBP-19 are significantly upregulated in SB. For OBP-C12, a role in the defense of T. castaneum has been already described in .
A study in T. castaneum conjectured an interaction between immune and developmental genes . Therefore, the overrepresentation analysis was also applied to a list of developmental genes taken from ; see Additional file 13. As can be observed in Figure 5, the category “developmental genes” is significantly overrepresented in the significantly upregulated genes of the Cro1 samples BttP, 6 h after infection, PC, 6 h and 18 h after exposure. Furthermore, these developmental genes are enriched in the significantly downregulated genes of the Cro1 sample BttO, 18 h after infection. Remarkably, the developmental genes are not enriched in any of the SB samples, neither for up- nor for downregulated genes. Thus, up- and downregulation of developmental genes upon oral infection, pricking infection and pricking control seem to constitute an important response mechanism of the beetle population Cro1 while for SB, the involvement of developmental genes in response to these types of exposures appears to be negligible.
Immune gene categories show population-specific differential regulation
In this study, we compared for the first time host gene expression profiles for differing routes of infection with the same bacterium, B. thuringiensis. We found that the routes of infection, oral and systemic (by pricking the cuticle) induce strikingly different transcriptional profiles in the host. Secondly, by using two host populations, we could show that the host populations show specific differences in their transcriptional responses. Intriguingly, the factors infection-route and host population interact, such that population differences depend on the infection routes and were generally more pronounced for oral infection.
Genetic underpinnings of infection route-specificity
A number of pathogens are able to infect by different routes, but in-depth analyses of global host response profiles when pathogens enter the body in different ways are as of yet lacking. Insects are important model organisms for studying responses to pathogen entry. Making primarily use of the fruit fly Drosophila melanogaster, injection or pricking of the cuticle for application of bacteria has been the method of choice that has enabled the identification of the main regulatory pathways of innate immunity [18–22]. More recently, protocols for oral infection with natural bacterial pathogens are increasingly being employed for studying gut immunity [23–26]. However, prior studies have focused on either way of infection, and as far as we are aware, no study has as of yet compared in detail the host responses to both ways of infection for the same pathogen. Transcriptome analysis (RNA-seq) was here used as a powerful tool to compare responses of the red flour beetle T. castaneum to oral and pricking (i. e. systemic) infection with B. thuringiensis. Strong differences in gene expression profiles were observed for both infection routes, at both an early (6 h) and a later (18 h) time point after infection. Essentially, the overlap of genes between oral and pricking infection was strikingly low, which clearly supports the recent findings that evolutionary adaptation of D. melanogaster to a bacterial pathogen is contingent upon the infection route . It is particularly noteworthy that different host-pathogen systems were used in  and the present study, which points to the generality of the phenomenon of infection route-specificity. Moreover, it demonstrates that physiological responses (this study) and evolutionary trajectories  are connected, i. e. both types of studies are complementary.
In our study, vegetative cells were used for pricking, while spore toxin mixtures were used for oral infection. Although this might account for some differences in host gene expression, the specific processes following the two routes of infection are rather likely responsible for the strong differences. Opposite to systemic infection, B. thuringiensis infecting via the oral route has to use different strategies to enter the host. In contrast to the situation in the hemolymph, in the gut, B. thuringiensis has to prevent ejection by stopping or reducing peristalsis, and it has to breach the peritrophic membrane and gut epithelium. The different physiological consequences of oral and pricking infection are clearly represented in the strongly diverging GO term (Figure 4) and enrichment analyses (Figures 5 and 6). Pricking infection entails a strong signal for stress and cytochrome P450 genes, supporting a previous study . Differences in infection routes are particularly obvious for immune genes. Pricking infection provides the expected immune gene activation, resulting in induction of Toll and IMD pathways (Figures 6 and 7) that lead to strong AMP expression. Moreover, further execution mechanisms such as ROS and cellular responses are activated. Similar trends were also observed in a recent RNA-seq study in which Tenebrio molitor beetles were injected with heat-killed Staphylococcus aureus bacteria . By contrast, oral infection provides a more complex picture, where the two beetle populations showed surprisingly contrasting gene activation patterns, as will be outlined in more detail below.
Host populations differ in expression profiles upon oral infection
We used two host populations for our comparison of infection-route specific gene expression profiles: a widely used laboratory strain (SB) and a newly collected strain (Cro1). The latter strain was previously shown to be more resistant to oral infection . A comparison of gene expression upon oral infection may thus help to identify genetic underpinnings of the difference in resistance. Accordingly, we found much stronger transcriptional responses to infection in the Cro1 strain. Most interestingly, while the populations reacted rather uniformly to pricking infection, expression patterns of certain groups of genes showed markedly differing responses to oral infection. Enrichment of cuticular protein genes showed such population-specific patterns for oral infection (up-regulated in Cro1, while down-regulated in SB; Figure 5). By contrast, pricking did not lead to population-specific patterns for this class of genes. The strong pattern for cuticular protein genes could be a consequence of infection-induced developmental asynchrony compared to naive beetle larvae, since molting entails cuticular turnover, as suggested by , for honeybees 36 h after oral infection with Paenibacillus larvae, the cause of American Foulbrood Disease. However, this might be somewhat less likely in our study, since differences were observed already at an early time-point as short as 6 h post infection (p. i.). Accordingly, we found enrichment for down-regulated developmental genes in Cro1, but not in SB, only at the later time point of 18 h post oral infection (Figure 5). It is possible that Cro1 is more resistant because it shows developmental adjustment when exposed to B. thuringiensis. This hypothesis is supported by our recent demonstration of oral priming for resistance in T. castaneum, which also entailed developmental retardation . By contrast, T. castaneum that were prick-primed with heat-killed bacteria were previously shown to speed up development and thereby potentially escape a perceived risk of infection . This effect differed among beetle lines, which is supported by the observation in the present study that developmental genes were enriched among up-regulated genes in the pricking control in the Cro1, but not the SB population. It could be that the Cro1 population has retained an adaptive developmental flexibility that has been lost in the SB strain that was maintained in the laboratory for many generations, often subjected to a rigid developmental timing and in absence of relevant pathogen pressure. These differences among beetle lines are indicative of genetic variance in the interaction between immunity and development and are clearly a fascinating field of further investigation. Epigenetic regulation might connect immunity and development, as recently demonstrated for a lepidopteran host, the greater wax moth Galleria mellonella. An alternative explanation for the relevance of genes with cuticular functions is the relevance of the peritrophic membrane for oral infection. This inner lining of the midgut contains chitin, and degradation of the peritrophic membrane is an important step during gut infection, since it enables bacteria to gain access to epithelial cells . Peritrophic membrane repair and induction of chitin-specific enzymes can be expected during gut infection, and accordingly, chitin deacetylase was found up-regulated upon oral infection. Chitin deacetylases are important for chitin turnover in T. castaneum. Moreover, in the context of infection, it might be potentially interesting that chitin deacetylation results in the formation of chitosan, a powerful antimicrobial agent [32, 33].
Based on oral infection of D. melanogaster with Erwinia carotovora, a study from 2009 suggests that gut homeostasis is maintained through a balance between cell damage due to the collateral effects of bacteria killing and the repair of the gut barrier, i. e. peritrophic membrane and epithelia . Moreover, this analysis suggests that the IMD pathway directly participates in the remodeling of this barrier, while no evidence for the involvement of the Toll pathway in immune response against this oral infection was found . Our findings support this view, however, the situation seems even more complex, because the two host populations showed different responses (Figure 7). SB shows a pattern that is mostly consistent with D. melanogaster infected with E. carotovora, i. e. an up-regulation of the IMD pathway and therefore several AMPs. However, SB beetles also showed up-regulation of some (but not all) components of the Toll pathway. By contrast, Cro1 does not show a clear up-regulation of the IMD pathway and even a down-regulation of some of the Toll receptors. This relates to a clearly distinct pattern of the two populations regarding effector functions involved in oral infection: SB basically shows an AMP-based response to oral infection, while Cro1 seems to rely on immunity provided by PO and ROS (Figure 6). Arguably, the latter type of response might be more efficient to clear an oral infection with B. thuringiensis, leading to increased resistance of the Cro1 population. ROS provide a very fast response to infection . As the infection outcome of B. thuringiensis in T. castaneum is dose dependent , fast reduction of parasite load at the beginning of infection may indeed improve host survival. However, ROS may also lead to oxidative damage of self tissue and its need for renewal , i. e. there might be a cost the host has to pay for protection against infection. Accordingly, an enrichment among up-regulated genes was observed for stress genes 6 h p. i., followed by Cytochrome P450 (CYP) genes 18 h p. i., only in the Cro1, but not the SB population (Figure 5). CYP are a diverse class of enzymes with versatile functions including metabolic detoxification and resistance [35, 36]. The observed up-regulation of stress and CYP genes may thus function to alleviate immunity-related self-damage. Taken together, these data show that host populations may markedly differ in their way they deal with bacterial oral infections and its collateral damage. The SB population has been in the laboratory since many generations, likely with reduced exposure to pathogens, which might have led to the evolutionary loss of efficient, but costly immune strategies. Such rapid evolution is possible, as also demonstrated by . The observed population differences in immunity to oral infection might also result from differences in the gut microbiota of the host populations. Gut microbiota have for example been shown to be involved in immune priming for resistance against Plasmodium in mosquitoes . In this study, immune priming was related to responses of hemocytes to the natural gut microbiota that enters the hemocoel upon penetration of the gut epithelium by Plasmodium ookinetes. Since B. thuringiensis oral infection also damages the gut epithelium, it is not unlikely that part of the response to oral infection could be a response to gut microbiota. Potential differences in gut microbiota of the two study populations might even explain the observed differences in host responses to infection. Interestingly, gut microbiota have been suggested to be relevant in B. thuringiensis infections in lepidopteran hosts [38–40]. These alternative or complementary explanations will be a fascinating field for further research.
Genes involved in oral infection with B. thuringiensis
A number of genes have previously been reported to be involved in the oral infection process of B. thuringiensis and are thus candidate resistance mechanisms. The Cry toxin of B. thuringiensis is the most intensively studied virulence factor. Of the reported Cry receptors [41, 42], several were found down-regulated upon oral exposure. Aminopeptidase N-like genes, E-cadherin (TcCad1) and sodium solute symporter protein (TcSSS) were down-regulated in both host populations at 6 h p. i. Although TcCad1 and TcSSS were found to bind the Cry3Ba more strongly than the Cry3Aa toxin of the beetle-infective B. thuringiensis strain tenebrionis that was used here , its down-regulation in our experiment may nevertheless suggest a role in infection with B. thuringiensis tenebrionis. Apolipophorin III (ApoLp-III) gene transcripts were elevated upon oral infection. This gene has been involved in host defense in several insect species [43, 44] and was recently shown to participate in the regulation of hemolymph PO activity in larvae infected with Cry3Ba producing strain . Further genes with a hypothesized role in defense against B. thuringiensis in T. castaneum were found up-regulated in our study [16, 42]. A variety of odorant binding proteins (OBPs) and chemosensory proteins (CSPs) were significantly up-regulated 6 h p. i. in both beetle populations. Generally, CSPs and OBPs are involved in olfaction and chemical communication in insects, but the functions of members of these large families of genes might be diverse [8, 45]. Elevated levels of CSP gene transcripts were previously reported after oral exposure to B. thuringiensis or its toxins in Tenebrio molitor and T. castaneum, and a role of OBP-C12 in T. castaneum defense against B. thuringiensis was also demonstrated . The beetle populations differed in the specific sets and the expression dynamics of these genes, some of which were also differentially expressed upon pricking infection. For example, whereas OBP genes generally showed up-regulation 6 h after oral infection, Cro1 showed down-regulation of OBP genes after pricking infection.
Finally, Osiris 18 and 19 were strongly up-regulated upon oral infection specifically in the Cro1 population while they were down-regulated in SB. Osiris comprises a family of highly conserved insect proteins of unknown function . Interestingly, Osiris genes were recently shown to be up-regulated in honeybees upon oral infection with Paenibacillus larvae although developmental asynchrony could not be excluded as an indirect cause of differential expression . Nevertheless, Osiris genes emerge as interesting targets for future studies of their potential role in defense against orally infecting pathogens.
Gene expression upon pricking infection contains a strong signal of wounding
Pricking infection consists of wounding and delivery of bacteria into the hemocoel. Therefore, transcriptional responses to pricking infection resemble wounding responses, but show additional components. This is clearly visible from differentially expressed gene numbers (Figure 3) and the GO term enrichment analyses (Figure 4), where pricking infection resembles wounding, while oral infection shows clearly distinct patterns. The strong overlap of differentially expressed gene numbers and GO term enrichment between bacterial pricking (BttP) and sterile wounding as a pricking control (PC) shows that the wounding response makes up for a large part of the observed effects of bacterial pricking. Many genes overrepresented for both BttP and PC treatment were identified as cuticle related in response to the physical damage of pricking and potentially also the induced developmental asynchrony of both treatments compared to naive beetles.
Also regarding stress related genes, Cytochrome P450, odorant binding and developmental genes, the bacterial pricking and sterile wounding treatments show rather similar enrichment (Figure 5). However, a more detailed comparison of immune pathways also revealed genes that are up-regulated specifically after bacterial pricking but not in the pricking control, such as Dif2, IMD and REL2 (Relish), which are central to the regulation of the two main immune-inducible pathways IMD and Toll; see Additional file 14. This might indicate that bacterial pricking (compared to wounding alone) leads to a stronger or more long-lasting activation of these pathways. While IMD is mostly involved in defense against Gram-negative bacteria, Toll is generally activated by Gram-positive bacteria and fungi (see  for a review on Drosophila immunity and  for T. castaneum). However, even though Bacillus is Gram-positive, it has DAP-type peptidoglycan, which is characteristic of Gram-negatives and activates the IMD pathway. Accordingly, we found both pathways activated after pricking infection with B. thuringiensis, and a large number of AMPs including Attacin1, Attacin2, Cecropin2, Cecropin3, Defensin1, Defensin3, Lysozyme2 and Lysozyme3 were strongly induced. In sum, our study clearly shows that sterile wounding is a necessary control to distinguish, which genes are involved in the wounding, and which are additionally relevant for combating the bacteria.
Our study shows that different routes of infection lead to strongly divergent gene expression profiles in the host, and that it is therefore important to include the aspect of infection routes into studies of host-pathogen interactions. Importantly, host populations differed in their responses in particular to oral infection. This suggests a higher degree of phenotypic and genetic variance for responses to oral infection, as compared to the more ‘hard-wired’ responses to pricking infection. It thus also corroborates a recent finding that evolutionary adaptation was faster for resistance to oral as compared to pricking infection . Host population differences emerge as an interesting field for future investigation, since a recently collected beetle population (Cro1) seems to react with more powerful immune reactions and to alleviate collateral damage, while a laboratory strain (SB) seems to have lost this flexibility. Resistance mechanisms showing strong diversity among populations are presumably evolving fast, and there is indeed some overlap of the population-specific resistance mechanisms identified in our study and the degree of adaptive evolution detected in immune system genes in Drosophila. Moreover, we found indications that the more outbred population may show a reduced degree of adaptive developmental plasticity. Genetic diversity in the interaction of immunity with development might thus be a relevant field for future studies to advance our understanding of resistance in an ecological and evolutionary context.
Two populations of T. castaneum were examined in this study. As a standard laboratory strain San Bernardino (SB) was used. Furthermore, a wild type strain, Croatia 1 (Cro1), collected in June 2010 in Croatia, was included . This strain was adapted to lab condition for more than 20 generations (18 months). Beetles were reared on flour (type 550) with 5% brewer yeast at 30°C with a 12/12 h light/dark cycle.
Bacillus thuringiensis bv. tenebrionis (Btt) spores were purchased from the Bacillus Genetic Stock Center (BGSC) and subcloned five times on LB-Agar before a log-phase culture was used for glycerin stocks that were stored at −80°C.
Approximately 500 two weeks old adults of each strain were used for egg production for 24 h, respectively. Eggs were sieved of the flour with 280 μm sieve and further cultivated under given standard conditions for another 10 days. Larvae were sieved of and individualized into 96 well plates containing flour and 5% yeast. After an additional time span of six days larvae of each strain were randomly assigned to the following treatments: Btt pricking (BttP), Btt oral (BttO), pricking control (PC), naïve control (NC). For analysis of differential RNA expression over time we choose 2 time points (6 h; 18 h) to sample total RNA. Each treatment was done with both beetle strains and 3 replicates with 32 larvae each (2 beetle strains × 4 treatments × 2 time points × 3 replicates = 48 samples). For BttP, bacteria were grown from a glycerin stock over night for 15 h in 50 ml standard LB-Media and centrifuged for 15 min, 5000×g at 4°C and washed with PBS twice. After last wash the pellet was resuspended in 2 ml PBS and counted in a Thoma counting chamber to adjust the concentration to 1 × 1010 per ml (LD20). For infection the larva was pricked dorsally between 1st and 2nd integument into the main vessel with a sterile dissecting needle (Ø10 μ m) that was dipped into either the bacteria solution (BttP) or into PBS (PC). After pricking, larvae were individualized in 96 well plates containing flour and 5% yeast until sampling. For BttO spore production, preparation of spore-containing diet and the infection protocol were done as previously described . Briefly, spores were produced and the suspension was centrifuged at 2880 × g at RT for 15 min, washed once, resuspended in PBS and counted using a Thoma counting chamber. Freshly produced spores were mixed with the beetle diet (0.15 g of flour with yeast/mL) in a concentration of 1 × 109 spores per ml (LD20) of diet. 40 μl of spore-containing liquid diet was pipetted into each well of a 96 well plate and dried overnight at 50°C. For infection, larvae were exposed to the spore-containing diet for three hours after which they were transferred to the spore-free diet until sampling. This was done since longer exposure does not contribute to further mortality to a higher extent . Control larvae were treated in the same way. Larvae were sampled 6 and 18 h after the initial exposure had started (zero time point being the time point when larva was placed on the disc), since the majority of larvae usually start feeding within minutes after being placed on the discs (personal observation). For each NC sample 16 larvae were kept under the same conditions as in BttP and PC and another 16 larvae were kept under the same conditions as in BttO without infection. For RNA isolation, 32 larvae of each treatment were pooled for one sample after 6 h and 18 h, respectively. Samples were shock frozen in liquid nitrogen and stored at −80°C. Total RNA from frozen beetles was isolated using mirVana TM miRNA Isolation Kit (Ambion) according to the instructions of the manufacturer; for an overview of the infection experiment see Figure 1.
Library preparation, sequencing and mapping
The libraries were created with the Illumina TruSeq RNA Sample Prep Kit v2. After cluster generation with the TruSeq PE Cluster Kit v3 (cBot –HS) the sequencing was performed with the TruSeq SBS Kit v3 –HS (200 cycles) on the Illumina HiSeq 2000. Before the mapping, the paired end reads of length 101 bp were preprocessed in the following way: (1) adapter sequences were removed from the data set using SeqPrep , (2) reads that did not pass the internal Illumina quality filter were eliminated and (3) due to a non-uniform distribution of nucleotides at the beginning of the reads – a problem naturally occurring caused by the random hexamer priming step in cDNA generation  – the first 13 base pairs of every read were trimmed using the FASTX-toolkit, version 0.0.13 . Afterwards, the trimmed and filtered paired-end reads were mapped against the T. castaneum reference genome, version 3.0, downloaded from Beetlebase [53–55], making use of the RNA-seq read mapper Tophat, version 2.0.8b , with Bowtie version 2.1.0 . Setting option -G, Tophat was also supplied with the T. castaneum genome annotation in gtf format, downloaded from the iBeetle webpage in June 2013 (“AUGUSTUS 2 prediction”) [58, 59]. iBeetle is an initiative of the group around Prof. Mario Stanke to enhance the current Beetlebase annotation of the T. castaneum genome by integrating data of various RNA sequencing experiments. In the “Beetle community”, this unofficial genome annotation is considered to be more accurate compared to the current official T. castaneum annotation from Beetlebase.
Transcript assembly, library normalization and differential expression analysis
For every sample, transcript assembly and quantification have been carried out using cufflinks from Cufflinks, version 2.1.1, with option -G which tells Cufflinks to use the supplied reference gft annotation file from iBeetle and with option --upper-quartile-norm [58–60]. This latter option entails an upper quartile (75th percentile) normalization within each library and improves sensitivity without loss of specificity for differential expression calls . Then the Cufflinks utility cuffmerge from Cufflinks was run to merge these assemblies into a comprehensive transcriptome. Significant changes in transcript expression between two different conditions were then detected employing the Cufflinks tool cuffdiff with option --upper-quartile-norm and the default q-value cutoff of 0.05. The utility cuffdiff performs this upper-quartile normalization across the whole set of samples and thus accounts for differences in library size and sequencing depth. To assess the distributions of upper-quartile normalization based FPKM (Fragments Per Kilobase of exon per Million fragments mapped) values across samples, a boxplot was produced using the R package cummeRbund ; see Additional file 15.
Principal component analyses
The principle component analyses were based on the normalized FPKM values calculated with cuffdiff and performed with the method pca in R (package “labdsv”) using the covariance matrix . The scores of the first two dimensions were plotted with the R function plot.
Gene annotation and functional analysis
Genes of interest were annotated in two different ways. First, iBeetle genes were matched against the official Beetlebase gene identifiers by using a best reciprocal Blast hits based association table kindly provided by the iBeetle consortium [58, 59]. Using the official Beetlebase gene identifiers, gene descriptions, Gene Ontology (GO) terms  and InterPro attributes were downloaded from the EnsemblMetazoa database, release 17 . In order to further improve the annotation of iBeetle genes, especially for those without a best reciprocal Blast hit among the Beetlebase gene identifiers, peptide sequences for genes of interest were downloaded from the iBeetle webpage and were blasted against the nr protein database using blastp from Blast2GO, version 2.6.6, and were then mapped and annotated with the Blast2GO default parameters [58, 59, 66]. Additionally, the Blast2GO InterProScan was applied and the InterProScan GOs were merged to the annotation .
Overrepresentation analysis of GO terms and TermLogos
In order to detect GO terms assigned by Blast2GO that are overrespresented in a gene subsample of interest against all T. castaneum genes as background, the one-tailed Fisher’s Exact Test implemented in Blast2GO with a Benjamini Hochberg corrected p-value cutoff of 0.05 was applied . Similarly to , using a scaling factor of S=| log10(p-value)| and a color scheme of length three to differentiate between the three different ontologies “biological process”, “cellular component” and “molecular function”, a TermLogo was generated making use of the web server Wordle T M . The resulting word cloud thus visualizes GO terms that are enriched in a subsample of interest against all T. castaneum genes as background such that the color within the TermLogo represents one of the three ontologies and the font size the significance of every GO term according to the overrepresentation analysis: the larger the font, the smaller the corresponding p-value .
Overrepresentation analysis of gene categories
Based on a list of 388 T. castaneum immunity-related genes published by  (of which 305 genes could be matched to genes predicted by iBeetle; of those genes, 303 genes could be found among assembled transcripts in our RNA-seq data set) separated into the categories “recognition”, “extracellular signaling”, “intracellular signaling” (with subcategories “Toll Pathway”, “IMD pathway”, “JNK pathway” and “JAK-STAT pathway”) and “execution” (with subcategories “PO”, “ROS”, “AMP”, “cellular responses”), we tested whether any of these (sub-) categories is overrepresented in the list of either significantly up- or downregulated genes (q-value cutoff: 0.05) for every pairwise Cufflinks differential expression result using Fisher’s Exact Test. For example, we tested whether there is a significant overrepresentation of “AMP” genes in the list of significantly upregulated genes resulting from the comparison of Cro1:BttO:6h vs. Cro1:NC:6h. The resulting p-values were corrected for multiple testing according to the Benjamini Hochberg procedure. Afterwards, the adjusted p-values were visualized using the function heatmap.2 from the R library “gplots” . Apart from the list of immune genes by , this enrichment analysis was also applied to additional genes sets. We defined a set of cuticle genes based on the T. castaneum GO annotation downloaded from the EnsemblMetazoa database, release 17 by extracting only those genes with the GO term “structural constituent of cuticle”. The resulting list of cuticle genes contains 120 genes of which 111 could be linked to iBeetle gene predictions. A further list of stress-responsive genes involving heat shock genes was defined by taking all T. castaneum genes that contain the string “stress” or “heat shock” in at least one GO term downloaded from the EnsemblMetazoa database. Of this list of 58 stress genes, 51 could be matched to iBeetle predictions. Furthermore, a set of P450 genes was defined by extracting all T. castaneum genes with the term “P450” in the corresponding EnsemblMetazoa gene’s description (120 genes of which 89 have an iBeetle equivalent). Similarly, the list of odorant binding genes is based on extracting all T. castaneum genes with the GO term “odorant binding” resulting in 262 genes of which 103 could be linked to iBeetle predictions. A list of developmental genes was defined based on Table S11 from  containing 442 “selected developmental genes” of which 397 genes could be both matched to genes predicted by iBeetle and found among assembled transcripts in our RNA-seq dataset. After applying Fisher’s Exact Test to these gene lists, the resulting p-values were corrected for multiple testing and visualized in heatmaps as stated above.
Availability of supporting data
The data sets supporting the results of this article are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA), accession number SRP033773.
This work was supported by grants BO 2544/7-1, KU 1929/4-1, RO 2994/3-2, SCHU 1415/9-2 within the DFG (German Research Foundation, http://dfg.de/) priority programme 1399 “Host parasite Coevolution”. Furthermore, we acknowledge support by the DFG and Open Access Publication Fund of University of Münster. We thank the iBeetle consortium, in particular Prof. Gregor Bucher and the group of Prof. Mario Stanke, and the following collaborators for contribution of unpublished RNA and EST data: Prof. Ernst Wimmer, Prof. Reinhard Schröder, Dr. Michael Schoppmeier, Prof. Andreas Vicinskas, Dr. Boran Altincicek, Prof. Herbert Jäckle, Dr. Ulrike Löhr, Dr. Ho-Ryun Chung, Prof. Sue Brown, Prof. Dick Beeman, Prof. Yoonseong Park, Dr. Marce Lorenzen, Prof. Markus Friedrich. Additionally, we thank Dr. Antoine Branca, Dr. Sonja Grath, Dr. Jennifer Greenwood and Dr. Mahesh Panchal for helpful discussions. Erich Bornberg-Bauer’s ORCID number is http://orcid.org/0000-0002-1826-3576.
- Vizoso DB, Ebert D: Phenotypic plasticity of host-parasite interactions in response to the route of infection. J Evol Biol. 2005, 18 (4): 911-921.PubMedView ArticleGoogle Scholar
- Leggett HC, Cornwallis CK, West SA: Mechanisms of pathogenesis, infective dose and virulence in human parasites. PLoS Pathog. 2012, 8 (2): 1002512-View ArticleGoogle Scholar
- Martins NE, Faria VG, Teixeira L, Magalhaes S, Sucena E: Host adaptation is contingent upon the infection route taken by pathogens. PLoS Pathog. 2013, 9 (9): 1003601-View ArticleGoogle Scholar
- Milutinović B, Stolpe C, Armitage SA, Kurtz J, Peuß R: The red flour beetle as a model for bacterial oral infections. PLoS ONE. 2013, 8 (5): 64638-View ArticleGoogle Scholar
- Roth O, Sadd BM, Schmid-Hempel P, Kurtz J: Strain-specific priming of resistance in the red flour beetle, Tribolium castaneum. Proc Biol Sci. 1654, 276: 145-151.View ArticleGoogle Scholar
- Wang Z, Gerstein M, Snyder M: RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009, 10 (1): 57-63.PubMed CentralPubMedView ArticleGoogle Scholar
- Klingler M: Tribolium. Curr Biol. 2004, 14 (16): 639-640.View ArticleGoogle Scholar
- Richards S, Gibbs RA, Weinstock GM, Brown SJ, Denell R, Beeman RW, Gibbs R, Beeman RW, Brown SJ, Bucher G, Friedrich M, Grimmelikhuijzen CJ, Klingler M, Lorenzen M, Richards S, Roth S, Schröder R, Tautz D, Zdobnov EM, Muzny D, Gibbs RA, Weinstock GM, Attaway T, Bell S, Buhay CJ, Chandrabose MN, Chavez D, Clerk-Blankenburg KP, Cree A, Dao M, et al: The genome of the model beetle and pest Tribolium castaneum. Nature. 2008, 452 (7190): 949-955.PubMedView ArticleGoogle Scholar
- Schnepf E, Crickmore N, Van Rie J, Lereclus D, Baum J, Feitelson J, Zeigler DR, Dean DH: Bacillus thuringiensis and its pesticidal crystal proteins. Microbiol Mol Biol Rev. 1998, 62 (3): 775-806.PubMed CentralPubMedGoogle Scholar
- van Frankenhuyzen K: Insecticidal activity of Bacillus thuringiensis crystal proteins. J Invertebr Pathol. 2009, 101 (1): 1-16.PubMedView ArticleGoogle Scholar
- Roth O, Joop G, Eggert H, Hilbert J, Daniel J, Schmid-Hempel P, Kurtz J: Paternally derived immune priming for offspring in the red flour beetle, Tribolium castaneum. J Anim Ecol. 2010, 79 (2): 403-413.PubMedView ArticleGoogle Scholar
- Moore AD, Bornberg-Bauer E: The dynamics and evolutionary potential of domain loss and emergence. Mol Biol Evol. 2012, 29 (2): 787-796.PubMed CentralPubMedView ArticleGoogle Scholar
- Johnston PR, Makarova O, Rolff J: Inducible defenses stay up late: temporal patterns of immune gene expression in Tenebrio molitor. G3 (Bethesda). 2013, doi:g3.113.008516Google Scholar
- Zou Z, Evans JD, Lu Z, Zhao P, Williams M, Sumathipala N, Hetru C, Hultmark D, Jiang H: Comparative genomic analysis of the Tribolium immune system. Genome Biol. 2007, 8 (8): 177-View ArticleGoogle Scholar
- Altincicek B, Elashry A, Guz N, Grundler FM, Vilcinskas A, Dehne HW: Next generation sequencing based transcriptome analysis of septic-injury responsive genes in the beetle Tribolium castaneum. PLoS ONE. 2013, 8 (1): 52004-View ArticleGoogle Scholar
- Contreras E, Rausell C, Real MD: Proteome response of Tribolium castaneum larvae to Bacillus thuringiensis toxin producing strains. PLoS ONE. 2013, 8 (1): 55330-View ArticleGoogle Scholar
- Roth O, Kurtz J: The stimulation of immune defence accelerates development in the red flour beetle (Tribolium castaneum). J Evol Biol. 2008, 21 (6): 1703-1710.PubMedView ArticleGoogle Scholar
- Boman HG, Nilsson I, Rasmuson B: Inducible antibacterial defence system in Drosophila. Nature. 1972, 237 (5352): 232-235.PubMedView ArticleGoogle Scholar
- Hoffmann JA, Reichhart JM: Drosophila innate immunity: an evolutionary perspective. Nat Immunol. 2002, 3 (2): 121-126.PubMedView ArticleGoogle Scholar
- Rolff J, Siva-Jothy MT: Invertebrate ecological immunology. Science. 472, 301 (5632):Google Scholar
- Lemaitre B, Hoffmann J: The host defense of Drosophila melanogaster. Annu Rev Immunol. 2007, 25: 697-743.PubMedView ArticleGoogle Scholar
- Chambers MC, Schneider DS: Pioneering immunology: insect style. Curr Opin Immunol. 2012, 24 (1): 10-14.PubMedView ArticleGoogle Scholar
- Nehme NT, Liegeois S, Kele B, Giammarinaro P, Pradel E, Hoffmann JA, Ewbank JJ, Ferrandon D: A model of bacterial intestinal infections in Drosophila melanogaster. PLoS Pathog. 2007, 3 (11): 173-View ArticleGoogle Scholar
- Vlisidou I, Dowling AJ, Evans IR, Waterfield N, ffrench-Constant RH, Wood W: Drosophila embryos as model systems for monitoring bacterial infection in real time. PLoS Pathog. 2009, 5 (7): 1000518-View ArticleGoogle Scholar
- Buchon N, Broderick NA, Poidevin M, Pradervand S, Lemaitre B: Drosophila intestinal response to bacterial infection activation of host defense and stem cell proliferation. Cell Host Microbe. 2009, 5 (2): 200-211.PubMedView ArticleGoogle Scholar
- Chakrabarti S, Liehl P, Buchon N, Lemaitre B: Infection-induced host translational blockage inhibits immune responses and epithelial renewal in the Drosophila gut. Cell Host Microbe. 2012, 12 (1): 60-70.PubMedView ArticleGoogle Scholar
- Cornman RS, Lopez D, Evans JD: Transcriptional response of honey bee larvae infected with the bacterial pathogen Paenibacillus larvae. PLoS ONE. 2013, 8 (6): 65424-View ArticleGoogle Scholar
- Fritzlar S, Kurtz J, Milutinović B: Increased survival in the red flour beetle after oral priming with bacteria-conditioned media. J Innate Immun. 2014, 6 (3): 306-314.PubMedGoogle Scholar
- Mukherjee K, Fischer R, Vilcinskas A: Histone acetylation mediates epigenetic regulation of transcriptional reprogramming in insects during metamorphosis, wounding and infection. Front Zool. 2012, 9 (1): 25-PubMed CentralPubMedView ArticleGoogle Scholar
- Lehane MJ: Peritrophic matrix structure and function. Annu Rev Entomol. 1997, 42: 525-550.PubMedView ArticleGoogle Scholar
- Arakane Y, Dixit R, Begum K, Park Y, Specht CA, Merzendorfer H, Kramer KJ, Muthukrishnan S, Beeman RW: Analysis of functions of the chitin deacetylase gene family in Tribolium castaneum. Insect Biochem Mol Biol. 2009, 39 (5–6): 355-365.PubMedView ArticleGoogle Scholar
- Li B, Su T, Chen X, Liu B, Zhu B, Fang Y, Qiu W, Xie G: Effect of chitosan solution on the bacterial septicemia disease of Bombyx mori (Lepidoptera: Bombycidae) caused by Serratia marcescens. Appl Entomol Zool. 2010, 45 (1): 145-152.View ArticleGoogle Scholar
- Kong M, Chen XG, Xing K, Park HJ: Antimicrobial properties of chitosan and mode of action: a state of the art review. Int J Food Microbiol. 2010, 144 (1): 51-63.PubMedView ArticleGoogle Scholar
- Cerenius L, Soderhall K: The prophenoloxidase-activating system in invertebrates. Immunol Rev. 2004, 198: 116-126.PubMedView ArticleGoogle Scholar
- Feyereisen R: Insect P450 enzymes. Annu Rev Entomol. 1999, 44: 507-533.PubMedView ArticleGoogle Scholar
- Scott JG, Wen Z: Cytochromes P450 of insects: the tip of the iceberg. Pest Manag Sci. 2001, 57 (10): 958-967.PubMedView ArticleGoogle Scholar
- Rodrigues J, Brayner FA, Alves LC, Dixit R, Barillas-Mury C: Hemocyte differentiation mediates innate immune memory in Anopheles gambiae mosquitoes. Science. 2010, 329 (5997): 1353-1355.PubMed CentralPubMedView ArticleGoogle Scholar
- Broderick NA, Raffa KF, Handelsman J: Midgut bacteria required for Bacillus thuringiensis insecticidal activity. Proc Natl Acad Sci U S A. 2006, 103 (41): 15196-15199.PubMed CentralPubMedView ArticleGoogle Scholar
- Broderick NA, Robinson CJ, McMahon MD, Holt J, Handelsman J, Raffa KF: Contributions of gut bacteria to Bacillus thuringiensis-induced mortality vary across a range of Lepidoptera. BMC Biol. 2009, 7: 11-PubMed CentralPubMedView ArticleGoogle Scholar
- Johnston PR, Crickmore N: Gut bacteria are not required for the insecticidal activity of Bacillus thuringiensis toward the tobacco hornworm, Manduca sexta. Appl Environ Microbiol. 2009, 75 (15): 5094-5099.PubMed CentralPubMedView ArticleGoogle Scholar
- Pigott CR, Ellar DJ: Role of receptors in Bacillus thuringiensis crystal toxin activity. Microbiol Mol Biol Rev. 2007, 71 (2): 255-281.PubMed CentralPubMedView ArticleGoogle Scholar
- Contreras E, Rausell C, Real MD: Tribolium castaneum Apolipophorin-III acts as an immune response protein against Bacillus thuringiensis Cry3Ba toxic activity. J Invertebr Pathol. 2013, 113 (3): 209-213.PubMedView ArticleGoogle Scholar
- Whitten MM, Tew IF, Lee BL, Ratcliffe NA: A novel role for an insect apolipoprotein (apolipophorin III) in beta-1,3-glucan pattern recognition and cellular encapsulation reactions. J Immunol. 2004, 172 (4): 2177-2185.PubMedView ArticleGoogle Scholar
- Gupta L, Noh JY, Jo YH, Oh SH, Kumar S, Noh MY, Lee YS, Cha SJ, Seo SJ, Kim I, Han YS, Barillas-Mury C: Apolipophorin-III mediates antiplasmodial epithelial responses in Anopheles gambiae (G3) mosquitoes. PLoS ONE. 2010, 5 (11): 15410-View ArticleGoogle Scholar
- Pelosi P, Calvello M, Ban L: Diversity of odorant-binding proteins and chemosensory proteins in insects. Chem Senses. 2005, 30 Suppl 1: 291-292.View ArticleGoogle Scholar
- Oppert B, Dowd SE, Bouffard P, Li L, Conesa A, Lorenzen MD, Toutges M, Marshall J, Huestis DL, Fabrick J, Oppert C, Jurat-Fuentes JL: Transcriptome profiling of the intoxication response of Tenebrio molitor larvae to Bacillus thuringiensis Cry3Aa protoxin. PLoS ONE. 2012, 7 (4): 34624-View ArticleGoogle Scholar
- Shah N, Dorer DR, Moriyama EN, Christensen AC: Evolution of a large, conserved, and syntenic gene family in insects. G3 (Bethesda). 2012, 2 (2): 313-319.View ArticleGoogle Scholar
- Yokoi K, Koyama H, Ito W, Minakuchi C, Tanaka T, Miura K: Involvement of NF-κB transcription factors in antimicrobial peptide gene induction in the red flour beetle, Tribolium castaneum. Dev Comp Immunol. 2012, 38 (2): 342-351.PubMedView ArticleGoogle Scholar
- Obbard DJ, Welch JJ, Kim KW, Jiggins FM: Quantifying adaptive evolution in the Drosophila immune system. PLoS Genet. 2009, 5 (10): 1000698-View ArticleGoogle Scholar
- St. John J: SeqPrep. [https://github.com/jstjohn/SeqPrep]
- Hansen KD, Brenner SE, Dudoit S: Biases in Illumina transcriptome sequencing caused by random hexamer priming. Nucleic Acids Res. 2010, 38 (12): 131-View ArticleGoogle Scholar
- Hannon laboratory: Fastx Toolkit. [http://hannonlab.cshl.edu/fastx_toolkit/]
- Kim HS, Murphy T, Xia J, Caragea D, Park Y, Beeman RW, Lorenzen MD, Butcher S, Manak JR, Brown SJ: BeetleBase in 2010: revisions to provide comprehensive genomic information for Tribolium castaneum. Nucleic Acids Res. 2010, 38 (Database issue): 437-442.View ArticleGoogle Scholar
- Wang L, Wang S, Li Y, Paradesi MS, Brown SJ: BeetleBase: the model organism database for Tribolium castaneum. Nucleic Acids Res. 2007, 35 (Database issue): 476-479.View ArticleGoogle Scholar
- K-State Bioinformatics Center: BeetleBase Website. [http://www.Beetlebase.org]
- Trapnell C, Pachter L, Salzberg SL: TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009, 25 (9): 1105-1111.PubMed CentralPubMedView ArticleGoogle Scholar
- Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10 (3): 25-View ArticleGoogle Scholar
- iBeetle consortium: iBeetle Website. [http://bioinf.uni-greifswald.de/tcas/index.html]
- Stanke M, Diekhans M, Baertsch R, Haussler D: Using native and syntenically mapped cDNA alignments to improve de novo gene finding. Bioinformatics. 2008, 24 (5): 637-644.PubMedView ArticleGoogle Scholar
- Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L: Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010, 28 (5): 511-515.PubMed CentralPubMedView ArticleGoogle Scholar
- Bullard JH, Purdom E, Hansen KD, Dudoit S: Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinformatics. 2010, 11: 94-PubMed CentralPubMedView ArticleGoogle Scholar
- Goff L, Trapnell C, Kelley D: cummeRbund: analysis, exploration, manipulation, and visualization of Cufflinks high-throughput sequencing data. R package version 2.0.0. [http://compbio.mit.edu/cummeRbund/]
- Roberts DW: labdsv: Ordination and Multivariate Analysis for Ecology. R package version 1.6-1. [http://CRAN.R-project.org/package=labdsv]
- Gaudet P, Chisholm R, Berardini T, Dimmer E, Engel SR, Fey P, Hill DP, Howe D, Hu JC, Huntley R, Khodiyar VK, Kishore R, Li D, Lovering RC, McCarthy F, Ni L, Petri V, Siegele DA, Tweedie S, Van Auken K, Wood V, Basu S, Carbon S, Dolan M, Mungall CJ, Dolinski K, Thomas P, Ashburner M, Blake JA, Cherry JM, et al: The Gene Ontology’s Reference Genome Project: a unified framework for functional annotation across species. PLoS Comput Biol. 2009, 5 (7): 1000431-View ArticleGoogle Scholar
- Kinsella RJ, Kahari A, Haider S, Zamora J, Proctor G, Spudich G, Almeida-King J, Staines D, Derwent P, Kerhornou A, Kersey P, Flicek P: Ensembl BioMarts: a hub for data retrieval across taxonomic space. Database (Oxford). 2011, 2011: 030-View ArticleGoogle Scholar
- Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M: Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005, 21 (18): 3674-3676.PubMedView ArticleGoogle Scholar
- Quevillon E, Silventoinen V, Pillai S, Harte N, Mulder N, Apweiler R, Lopez R: InterProScan: protein domains identifier. Nucleic Acids Res. 2005, 33 (Web Server issue): 116-120.View ArticleGoogle Scholar
- Warnes GR, Bolker B, Bonebakker L, Gentleman R, Huber W, Liaw A, Lumley T, Maechler M, Magnusson A, Moeller S, Venables B, Schwartz M: gplots: Various R Programming Tools for Plotting Data. R Package Version 2.11.0. [http://CRAN.R-project.org/package=gplots]
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.