Next Generation Sequencing Based Transcriptome Analysis of Septic-Injury Responsive Genes in the Beetle Tribolium castaneum

Beetles (Coleoptera) are the most diverse animal group on earth and interact with numerous symbiotic or pathogenic microbes in their environments. The red flour beetle Tribolium castaneum is a genetically tractable model beetle species and its whole genome sequence has recently been determined. To advance our understanding of the molecular basis of beetle immunity here we analyzed the whole transcriptome of T. castaneum by high-throughput next generation sequencing technology. Here, we demonstrate that the Illumina/Solexa sequencing approach of cDNA samples from T. castaneum including over 9.7 million reads with 72 base pairs (bp) length (approximately 700 million bp sequence information with about 30× transcriptome coverage) confirms the expression of most predicted genes and enabled subsequent qualitative and quantitative transcriptome analysis. This approach recapitulates our recent quantitative real-time PCR studies of immune-challenged and naïve T. castaneum beetles, validating our approach. Furthermore, this sequencing analysis resulted in the identification of 73 differentially expressed genes upon immune-challenge with statistical significance by comparing expression data to calculated values derived by fitting to generalized linear models. We identified up regulation of diverse immune-related genes (e.g. Toll receptor, serine proteinases, DOPA decarboxylase and thaumatin) and of numerous genes encoding proteins with yet unknown functions. Of note, septic-injury resulted also in the elevated expression of genes encoding heat-shock proteins or cytochrome P450s supporting the view that there is crosstalk between immune and stress responses in T. castaneum. The present study provides a first comprehensive overview of septic-injury responsive genes in T. castaneum beetles. Identified genes advance our understanding of T. castaneum specific gene expression alteration upon immune-challenge in particular and may help to understand beetle immunity in general.


Introduction
Parasites reduce the fitness of their hosts and therefore numerous host mechanisms have evolved to limit infectious diseases.In animals, the risk of an infection is reduced by physical and chemical barriers, by behavioral defense reactions such as avoidance or hygiene [1], and by the complex and highly evolved immune defense system.In vertebrates, the immune system is composed of the adaptive immunity including specific T-cell receptors and B-cell-derived antibodies and the evolutionarily more ancient innate immunity [2,3].Of note, vertebrate innate immunity shows many parallels to the invertebrate immunity.Insects, e.g.Drosophila melanogaster, have widely been used to elucidate invertebrate immune reactions.These reactions include entrapment of invading pathogens in clots, phagocytosis by immune-competent cells (hemocytes), and the induced production of antimicrobial peptides as well as reactive oxygen species, both underlying the induced expression of a wide array of immunerelated genes [4][5][6][7][8][9].
The recent determination of the Tribolium castaneum genome sequence [10] enabled the identification of numerous immunerelated genes by both homology-based [11] and experimental approaches [12].These studies provided first important insights into the T. castaneum immunity; however, our understanding of Tribolium immune responses is still fragmentary.The expression levels of only a limited number of Tribolium genes have been determined upon immune-challenge [11,12].To gain deeper insights into Tribolium immune responses, here, we investigated the whole transcriptome of naı ¨ve and immune-challenged beetles by Illumina/Solexa next generation sequencing.To induce strong immune responses in T. castaneum we used a commercially available crude lipopolysaccharide (LPS) preparation derived from Escherichia coli, which has widely been used as an elicitor of immune responses in numerous vertebrates and invertebrate species [12][13][14][15][16].
The present sequencing approach resulted in the identification of the transcriptome of T. castaneum and the identification of 70 genes with significantly elevated and 3 genes with reduced mRNA levels upon septic injury as determined by fitting the expression data with generalized linear models.

Biological samples for transcriptional analysis
The Tribolium stock that we used in this study was the T. castaneum wild-type strain San Bernardino.In contrast to the genome-sequenced GA-2 T. castaenum strain, the strain San Bernardino is ''wild-type'' since no consecutive generations of virgin single-pair, full-sib inbreeding were performed for 20 generations to obtain near-homozygous inbred condition needed for proper genome-sequencing [10].Beetles were maintained on whole-grain flour with 5% yeast powder at 31uC in darkness.For the experimental treatments, we have first randomly selected 40 young adult beetles (1-2 weeks after final ecdysis), which were subsequently divided by chance into two groups.LPS-challenge of 20 beetles was performed by ventrolaterally pricking of the imagoes abdomen using a dissecting needle dipped in an aqueous solution of 10 mg/ml lipopolysaccharide (LPS, purified Escherichia coli endotoxin 0111:B4, Cat.No.: L2630, Sigma, Taufkirchen, Germany), as described [12].At eight hours post LPS-challenge treated beetles and a biologically independent sample of 20 unstabbed, but similar handled beetles (control) were frozen in liquid nitrogen.We extracted total RNA from frozen beetles using the TriReagent isolation reagent (Molecular Research Centre, Cincinnati, OH, USA) according to the instructions of the manufacturer and synthesized cDNA samples using oligio-d(T) primers with the SMART PCR cDNA Synthesis Kit (Clontech, Mountain View, CA, USA) as previously described [12].Sequencing was done by the GATC GmbH (Konstanz, Germany) sequencing company on an Illumina GA2 sequencer.

Data analysis and bioinformatics
We have deposited the short read sequencing data with the following SRA accession numbers at NCBI sequence database: SRX022010 (immune-challenged beetles) and SRX021963 (naı ¨ve beetles).Sequencing reads were mapped by the sequencing company with ELAND Illumina software using the first 32 bp with highest sequencing quality and score values over 30 indicating 99.9% accuracy [17] and allowing one mismatch to the reference sequence of the Tribolium genome sequencing [18].To calculate statistical differences of the expression levels of genes between treatment and control and thereby to identify immuneresponsive genes we utilized DESeq package [19] within Bioconductor [20] and R [21].DESeq was used to normalize the count data, calculate mean values, fold changes, size factors, variance and P values (raw and adjusted) of a test for differential gene expression based on generalized linear models using negative binomial distribution errors.

Identification of Single Nucleotide Polymorphisms (SNPs) and Deletion Insertion Polymorphisms (DIPs) and de novo assembly
Single Nucleotide Polymorphisms (SNPs) and Deletion Insertion Polymorphisms (DIPs) detection tools within the CLC genomic workbench (version 4.9) were used to determine sequence variants.First, all Illumina reads were prepared by trimming of ambiguous nucleotides (.2 N) and low quality bases (,0.05).First we mapped all reads against the Glean assembly transcripts.Then, the level of SNPs and DIPs quality and significance was assessed by adjusting the quality filter to select only SNPs and DIPs that exists in a window of at least 11 bases and does not score more than 2 gaps or mismatches.The quality of the central base of each window was set to be at least 20 and the surrounding bases at least 15.The significance filter was adjusted to ignore SNPs and DIPs that have a coverage less than 4 and variant level less than 35% of corresponding reads.De novo assembly has been performed with the CLC genomics workbench (version 4.9) with the de novo assembly algorithm for Illumina reads with default parameters settings (Min.similarity allowed = 0.8 at length fraction = 0.5, deletion and insertion cost = 3, and mismatch cost = 2).

Sequence annotation
Sequence homology searches of predicted reference gene sequences (gleans) and subsequent functional annotation by gene ontology terms (GO), InterPro terms (InterProScan, EBI), enzyme classification codes (EC), and metabolic pathways (KEGG, Kyoto Encyclopedia of Genes and Genomes) were determined using the BLAST2GO software suite v2.3.1 [22].Homology searches were performed remotely on the NCBI server through QBLAST: sequences were compared with the NCBI non-redundant (nr) protein database and matches with an E-value cut-off of 10 23 , with predicted polypeptides of a minimum length of 15 amino acids, were scored.Subsequently, GO classification, including enzyme classification codes and KEGG metabolic pathway annotations, were generated.For final annotation, InterPro searches on the InterProEBI web server were performed remotely by utilizing BLAST2GO.

Mapping Illumina sequencing reads to predicted gene models of T. castaneum
To gain insights into Tribolium immune responses, we investigated the whole transcriptome of naı ¨ve and immune-challenged beetles by Illumina/Solexa next generation sequencing.This sequencing approach resulted in over 9.7 million cDNA reads with over 700 million bp sequence information and estimated 306 transcriptome coverage.About 3.8 and 4.0 million reads of Illumina sequencing of control and LPS-challenged animals, respectively, were mapped to predicted gene models of T. castaneum, which were built on the 3.0 genome assembly [10] (Table 1).We found that 11,679 predicted genes were expressed in both naı ¨ve and LPS-challenged adult Tribolium beetles.Additional sequences corresponding to the expression of further 642 and 739 predicted genes in naı ¨ve and LPS-challenged beetles, respectively, were also observed.In total, this approach resulted in the expression validation of 13,060 genes, representing almost 80% of the in total 16,422 predicted genes.
Evidence for the need of gene model curation and identification of single-nucleotide polymorphisms (SNPs) and DIPs (short deletion and insertion polymorphisms) About 14% of all sequencing reads could be assigned to published T. castaneum EST sequences or the genome sequence but not to predicted gene models indicating that several exons or genes might be miss-predicted in the current genome annotation.Therefore, we shared the present sequencing data with the beetleBase [23] and the iBeetle consortium [24], which are currently working on a next, more precise genome annotation.In addition, we identified over 155,000 positions of high quality single-nucleotide polymorphisms (SNPs) and 895 DIPs (short deletion and insertion polymorphisms) within the coding gene sequences between the T. castaneum strain San Bernardino used in the present analysis and the genome-sequenced strain Georgia GA-2 [10] (Table S1).This information might be helpful in future comparative studies investigating the potential impact of SNPs and DIPs on varying ecological traits of diverse T. castaneum strains.Furthermore, we performed a de novo assembly (Data S1), which might be helpful for future studies investigating e.g.alternative splicing events.
Interestingly, about 5% of all sequencing reads did not map to T. castaneum sequences but to sequences from other organisms such as the bacteria Escherichia coli, Bacillus subtilis, or Azotobacter vinelandi.These bacterial species may represent part of the beetle flora.

Validation of present Illumina sequencing approach by comparing estimated fold change expression values with recently reported values determined by qRT-PCR analysis
To determine differentially expressed genes between naive and LPS-challenged beetles we first checked whether sequencing samples were comparable.We counted the amount of reads aligned to predicted genes using only the first 32 bp of reads with highest sequencing quality and score values over 30 indicating 99.9% sequence accuracy [17] (Figure S1).In both treatments, we found that almost all genes were expressed at identical levels resulting in a significant linear correlation of the logarithmically transformed expression values (Figure 1).The regression analysis resulted in an adjusted R-squared value of 0.9073 (F, 1.143610 5 ; d.f., 11,677, P, ,2.2610 216 ).However, as expected, several potentially immune-responsive genes showed variance in their expression levels and we compared their expression rates with recently investigated immune-responsive genes [12].Validating our present approach, the expression values determined by our recent qRT-PCR analysis of the house-keeping gene a-tubulin as control and several antimicrobial peptides such as defensins and thaumatin as well as stress-responsive genes such as heat shock factors [12] were found to be comparable to the values determined by the present RNA-Seq approach (Table 2).We found that the values of both experiments were highly similar and correlated with statistical significance (Pearson correlation factor of 0.95 of logarithmically transformed values with a Holm's method adjusted P values = 0) (Figure 2).

Identification of significantly induced or repressed genes upon LPS-challenge in T. castaneum
To identify novel immune-responsive genes we calculated statistical differences of the expression levels between treatments utilizing DESeq package within Bioconductor and R.This powerful tool estimated the variance in our data and tested for differential gene expression [19].Since the two biological independent samples from control and treated beetles resulted in comparable expression values (F, 1.143610 5 ; d.f., 11,677, P, ,2.2610 216 ), we took the variance estimated from comparing their count rates across conditions as described in the DESeq manual [25].This analysis to identify differentially expressed genes is appropriate and will only cause the variance estimate to be too high, so that the test will err to the side of being too conservative [25].We further used pools of 20 individuals per sample to average across biological replicates of individuals.In sum, normalized count data were fitted with a generalized linear model (GLM) estimating a negative binomial distribution to the calculated mean values of the two biologically independent samples with each containing pooled cDNAs of 20 individual beetles.Then the P values were adjusted for multiple testing with the Benjamini-Hochberg procedure, which controls false discovery rate (FDR) (Table S2).Finally, we obtained the statistically significant up-regulation of 70 genes and down-regulation of 3 genes with a 5% FDR (Figure 3).
To assign the potential functions of identified genes we performed an annotation step with blast2go (Table S3) and summarized differentially expressed genes (Table 3).We observed the strongly induced expression of numerous genes including specific serine proteases, Toll receptor, or cathepsin L that are reportedly immune-responsive also in Drosophila flies [6,26].Moreover, we found several genes encoding proteins with leucine-rich-repeat domains potentially involved in immune signaling reactions in Tribolium, which have not been investigated yet.The leucine-rich repeat domain is a common structural motif for the molecular recognition of microbes, which is also present in the prominent Toll-like receptors, evolutionarily conserved receptors initiating signaling reactions in animal immunity [2].
Of note, we found that several genes encoding proteins with haemolymph juvenile hormone binding domains were significantly induced (e.g.Glean 09776 and 09775) while expression of a paralogous gene was significantly reduced (Glean 13657) upon immune-challenge.These homologues genes may regulate beetle developmental processes by influencing hormone levels.In agreement with this assumption, recent studies described significantly elevated metamorphosis rates [27] or accelerated aging rates [28] in immune-challenged beetles Two further significantly down-regulated genes encode proteins with one an esterasedomain and the other a heparin-binding domain both with unknown function.A deeper understanding of the molecular regulation of beetle development by immune responses would help to unravel potential ecological traits in Tribolium that might be traded-off with immune reactions probably similar as shown for other insects [29][30][31][32].Expression rates of immune-related genes upon LPSchallenge in T. castaneum The expression rates of numerous immune-related genes showed high induction levels, such as in the case of attacins and defensins (Table 4).However, due to the limitation of the present in-depth sequencing and calculation procedure, we observed statistical significance in immune-induced expression for only a limited number of immune-related genes (Table 3); short gene sequence and low expression rates of e.g.antimicrobial peptides in naı ¨ve animals resulted in a higher variance estimate and a lower confidence in the base mean estimates.Hence, only genes expressed both at medium or high rate and with at least more than 4 fold expression changes were identified by our approach (Figure 3).Particularly genes encoding antimicrobial peptides such as attacins or defensins are expressed at very low level in unchallenged beetles resulting in a high variance estimate in the present analysis resulting in much lower power of statistical analysis.To identify even more genes with significant expression difference a much higher coverage and more replicate determination per treatment with at least 3-fold deeper sequencing [33] would be needed.However, here we will compare tendencies of

Molecular pattern recognition proteins
Peptidoglycan recognition proteins (PGRPs) are evolutionarily conserved in animals and have been found to bind specifically to and to hydrolyze bacterial peptidoglycan.In addition, peptidoglycan-bound insect PGRPs activate the Toll and IMD signal transduction pathways as well as immune-related proteolytic cascades [4,5].Genome-wide gene expression profiling of the Drosophila immune-response implied that five PGRP genes including PGRP-SA, SC2, SB1, LB and SD are up regulated upon septic injury [6].Here, we found that 5 of 7 Tribolium genes encoding PGRP-SA, LA, LC, SB, and LB were up regulated in response to LPS injection whereas the expression rates of PGRP-LE and LD were not significantly influenced.In Drosophila, PGRP-LD is only expressed in hemocytes and its function is yet unknown Table 2. Comparison of RNA level estimation by our recent qRT-PCR analysis [12] and present transcriptome sequencing approach.whereas PGRP-LE is an intracellular receptor capable of binding bacteria in the cytoplasm [4].Gram negative bacteria binding proteins (GNBP) comprises a family of proteins, also known as b-1,3-glucan binding proteins (bGBP) and b-1,3-glucan recognition proteins (bGRP).The first b-1,3-glucan recognition protein was purified from the hemolymph of the silkworm Bombyx mori with a strong affinity for gram negative bacteria [34].This GNBP contained a region with significant sequence homology to the catalytic region of a group of bacterial b-1,3-glucanases.In Drosophila, three GNBP paralogs (GNBP1, GNBP2 and GNBP3) are known, which only GNBP3 (CG13422) is immune-responsive upon septic injury [4].GNBP1 is required for Toll activation in response to gram positive bacterial infection whereas GNBP3 has been reported to sense fungal infections [4].The biological function of Drosophila GNBP2 has yet not been determined.In Tribolium, we found up regulation of bGRP3 but not of bGRP1 and bGRP2 upon LPS-challenge, which resembles observations from Drosophila.
Thioester-containing proteins (TEPs) are a further group of bacteria-binding proteins, which function as both opsonins and protease inhibitors [4,5].In Drosophila, TEP II and TEP IV of in total 6 paralogous TEPs were found to be induced upon septic injury [6].In T. castaneum only 4 TEPs are traceable in the whole genome sequence and in this study, we found that mRNA levels of Tribolium TEP-B and C were increased but not of TEP-A and D upon immune-challenge.Interestingly, no clear orthologs can be assigned between dipteran and coleopteran TEPs, except for Tribolium TEP-A, which is orthologous to Drosophila TEP-VI [11].Finally, a putative TEP/complement-binding receptor-like protein (LpR2) was shown to be immune-inducible in Drosophila but failed to exhibit significant difference in gene expression rates in Tribolium in our approach.

Immunity related signaling
In insects, major immune-related signaling pathways include Toll, IMD, JAK-STAT, and JNK pathways [4].Mammalian Tolllike receptors are capable of directly binding danger (e.g.extracellular nucleic acids or uric acid) or pathogen-derived molecules (e.g.LPS) while Drosophila Toll is instead activated by a proteolytically activated cytokine-like molecule spaetzle 1.In D. melanogaster spaetzle 1 is immune-responsive [6] and five further paralogs have been described (spaetzle 2-6) with yet unknown functions.In Tribolium, 7 spaetzle paralogs exist with spaetzle 3, 4, 5, and 6 representing orthologs to respective Drosophila spaetzle isoforms.However, no Tribolium ortholog of Drosophila spaetzle 2 can be found and spaetzle 1, 2 and 7 in Tribolium form a clade together with the single, immune-responsive Drosophila spaetzle 1 [11].Here, we found that Tribolium spaetzle 2 and 7 are immuneinducible, whereas spaetzle 1 was immune-repressed.
Similarly to their potential ligands, also Toll-like receptors have experienced lineage-specific gene duplications in beetles as well as in flies or mosquitoes [11].4 Tribolium Toll-like receptors (1 to 4) of in total 10 paralogs were described to form a clade with the single, immune-responsive Drosophila Toll receptor of the in total 9 Drosophila Toll-like receptors [11].Here, we observed that these Tribolium Toll-like receptors 1 to 4 were induced in their gene expression upon LPS-challenge similar to the Drosophila Toll receptor upon septic injury [6].In addition, Tribolium Toll6 was over 2-fold down-regulated whereas Toll8 was slightly upregulated.Taken together, these results support the hypothesis that Tribolium has a more complex immune-related Toll signaling than Drosophila, since both a higher number of immune-responsive Tolllike receptors and of spaetzle ligands exist in Tribolium than in Drosophila.
Regarding further immune-related signaling pathways, we found 2 to 5 fold induced expression of several signaling proteins involved in IMD and JNK pathways such as IMD and D-Jun, respectively, which is in agreement to observations from Drosophila  [6].Also in agreement with observations in Drosophila [6], we found that expression rates of JAK-STAT pathway genes were not significantly influenced by LPS-challenge in Tribolium.

Antimicrobial peptides
As expected, we identified genes encoding antimicrobial peptides such as defensins, attacins and thaumatin among the systemically most septic injury inducible genes with up to 10 to 30 fold higher expression rates in LPS-challenged animals than in naive ones.This is in agreement with observations from diverse immune-challenged invertebrates [6,7,9,[35][36][37].

Stress response genes
Recently, we determined induced expression of genes in T. castaneum involved in detoxification and stress adaptation such as apolipoprotein D, cytochrome P450, gluthathione S-transferase, and a number of heat shock proteins [12].In line with these observations, here we found elevated gene expression rates of a number of stress and detoxification genes upon LPS-challenge including most notable HSPs, CytP450s (e.g.6BQ7, 345D2, 6BQ12, 6BK5), GST, ApoD, and ABC transporters (Table 3).This supports our recent hypothesis that interdependencies between immune and stress responses exist in T. castaneum [12,38].
It should be noted here that wounding itself can lead to gene expression alterations in insects triggered by e.g.cryptic, endogenous danger signals such as nucleic acids or collagen fragments [39,40].Moreover, the presently used LPS preparation is known to include bacterial nucleic acids and peptidoglycans, which may be responsible for the induction of e.g.PGRPs and PGRPcontrolled genes.Hence, in follow-up studies we propose to investigate transcriptomic immune responses from beetles with varying treatments such as feeding and stabbing with different elicitors and pathogens of diverse phylogenic origin and much more time points of samples derived from whole animals, specific organs, tissues or cell types.In addition, different T. castaneum genotypes, sexes or developmental stages are likely to vary in their immune investment and hence may show altered gene expression upon immune-challenge, particularly in the context of diverse environmental cues and stresses.

Conclusions
The beetle immune response underlies the differential expression of a wide array of different genes.Here we describe differential expression of numerous immune-related genes as well as several genes encoding proteins with leucine-rich-repeat domains, which might function as receptors in specific immune recognition and signaling reactions in beetles maybe in a similar way as leucine-rich-repeat domain containing receptors in ancient jawless vertebrates [41].While insect immune defense mechanisms had generally been assumed to be non-specific, diverse insects including the red flour beetle T. castaneum have recently been shown to respond quite specifically to some pathogens [42][43][44][45].Presently identified genes may help to elucidate the molecular basis of such specific reactions.This study is the first whole transcriptome analysis of the complex gene expression response in T. castaneum upon septic injury and provides numerous candidate genes that we can use as a starting point for further studies on beetle immunity.

Figure 1 .Figure 2 .
Figure 1.Gene expression in naive and immune-challenged beetles.All reads were aligned to predicted genes and are shown as log2 values derived from cDNA of naı ¨ve and LPS-challenged animals, respectively.The linear correlation is indicated by a red line (F-test, P, ,2.2610 216 ).doi:10.1371/journal.pone.0052004.g001

Figure 3 .
Figure 3. Significance plot.The log2 fold change value of each gene is shown against its base mean value.Differentially expressed genes with statistically significant difference at 5% FDR are indicated by red coloring.doi:10.1371/journal.pone.0052004.g003

Table 3 .
Transcripts with significant differential expression upon LPS-challenge in adult beetles.

Table 4 .
Expression levels of immune-related genes in adult beetles.