Identification of Immunity-Related Genes in Ostrinia furnacalis against Entomopathogenic Fungi by RNA-Seq Analysis

Background The Asian corn borer (Ostrinia furnacalis (Guenée)) is one of the most serious corn pests in Asia. Control of this pest with entomopathogenic fungus Beauveria bassiana has been proposed. However, the molecular mechanisms involved in the interactions between O. furnacalis and B. bassiana are unclear, especially under the conditions that the genomic information of O. furnacalis is currently unavailable. So we sequenced and characterized the transcriptome of O. furnacalis larvae infected by B. bassiana with special emphasis on immunity-related genes. Methodology/Principal Findings Illumina Hiseq2000 was used to sequence 4.64 and 4.72 Gb of the transcriptome from water-injected and B. bassiana-injected O. furnacalis larvae, respectively. De novo assembly generated 62,382 unigenes with mean length of 729 nt. All unigenes were searched against Nt, Nr, Swiss-Prot, COG, and KEGG databases for annotations using BLASTN or BLASTX algorithm with an E-value cut-off of 10−5. A total of 35,700 (57.2%) unigenes were annotated to at least one database. Pairwise comparisons resulted in 13,890 differentially expressed genes, with 5,843 up-regulated and 8,047 down-regulated. Based on sequence similarity to homologs known to participate in immune responses, we totally identified 190 potential immunity-related unigenes. They encode 45 pattern recognition proteins, 33 modulation proteins involved in the prophenoloxidase activation cascade, 46 signal transduction molecules, and 66 immune responsive effectors, respectively. The obtained transcriptome contains putative orthologs for nearly all components of the Toll, Imd, and JAK/STAT pathways. We randomly selected 24 immunity-related unigenes and investigated their expression profiles using quantitative RT-PCR assay. The results revealed variant expression patterns in response to the infection of B. bassiana. Conclusions/Significance This study provides the comprehensive sequence resource and expression profiles of the immunity-related genes of O. furnacalis. The obtained data gives an insight into better understanding the molecular mechanisms of innate immune processes in O. furnacalis larvae against B. bassiana.


Introduction
Innate immunity is an ancient and universal host defense mechanism found in both vertebrates and invertebrates. Insects lack adaptive immunity and mainly rely on innate immunity for defense against the invasion of pathogens or parasites [1,2]. Innate immunity in insects consists of a humoral response based on melanization and antimicrobial peptides (AMPs), and a cellular response including phagocytosis, encapsulation, or nodulation of pathogens [2][3][4][5][6]. Melanization, involving melanin synthesis, sequesters and kills invading microorganisms or parasites [7,8]. AMPs are directly cytotoxic to a wide range of microorganims, and are synthesized through the Toll and Imd (Immune deficiency) pathways [9][10][11][12].
Current understanding of insect innate immune mechanism is mainly from powerful genetic studies in Drosophila melanogaster [13,14], and biochemical studies in relatively large insects, such as the silkworm, Bombyx mori [15,16], the tobacco hornworm, Manduca sexta [17][18][19], and the beetle Tenebrio molitor [20,21]. Conceptually, the innate immune reactions in insects can be divided into three steps: (1) recognition of pathogen-associated molecular patterns; (2) signal modulation and transduction; (3) production or replenishment of immune-related molecules including effectors (for review see [22,23]). The recognition step is mediated by a group of proteins, known as pattern recognition proteins (PRPs) such as peptidoglycan recognition protein (PGRP), b-1,3-glucan recognition protein (bGRP)/gram-negative binding protein (GNBP), C-type lectin (CTL) ( [24] and references within). A proteolytic cascade composed of a series of serine proteases amplifies the initial recognition signal and triggers the activation of signaling pathways [22,25]. Ultimately, immune effectors including AMPs are induced in specific tissues, such as fat body (liver analogue) and hemocytes (insect blood cells). These effector molecules combat the invasive microorganisms in the hemolymph (insect blood), and play direct roles in insect innate immunity. So far, knowledge about the identification and involvement of immunity-related genes in various insects, especially in non-model insects, is still unclear and incomplete.
Recently, genome-wide analysis has helped to identify immunity-related genes and gene families in several insect species including 3 dipteran insects D. melanogaster [26], Anopheles gambiae [27], Aedes agypti [28], a hymenopteran insect Apis mellifera [29], a coleopteran insect Tribolium castaneum [30], and a lepidopteran insect B. mori [31]. However, the systematic analysis of immunityrelated genes is hampered in other insects whose genomic information is unavailable. The introduction of novel high through-put sequencing technologies provides insight into a comprehensive immune-gene repertoire of non-model insects. Second generation high-throughput DNA sequencing platforms such as Roche 454, SOLID and Illumina combined with the development of de novo assembly strategies have become powerful tools in transcriptomic studies for non-model organisms without a proper reference genome, and allow targeted identification of genes which are (differentially) expressed upon activation of immune responses [32,33]. This technology has been used, for example, to characterize the immunity-related genes in the beet armyworm Spodoptera exigua [34], the brown planthopper Nilaparvata lugens [35], and the wax moth Galleria mellonella [36] etc. It will greatly facilitate the future studies about analyzing the blueprint of host's genes under microbial challenge, especially for some important insect pests in which molecular information about the immunity-related genes is lacking.
The Asian corn borer, Ostrinia furnacalis (Guenée), is a major insect pest in Asia and causes serious damage on corn, sorghum, millet and cotton [37]. Control of this pest with chemical insecticides is currently hindered by the cryptic nature of larval behavior. Excessive use of chemical insecticides also leads to severe environmental pollution and insecticide residence. Therefore, entomopathogenic fungi become one of promising alternates for its control. The potential for suppression of O. furnacalis larvae by entomopathogenic fungus Beauveria bassiana has been proposed [38]. However, the molecular mechanisms involved in the interactions between O. furnacalis and B. bassiana are still largely unknown, especially under the conditions that the genomic information of O. furnacalis is absent currently. This will greatly restrict the further development and wider adoption of entomopathogenic fungi as control agents. The first step to resolve this problem could be comprehensive identification and characterization of immunity-related genes involved in the response of O. furnacalis larvae against B. bassiana. In this study, we used the Illumina sequencing and de novo assembly to explore the O. furnacalis immune response stimulated by B. bassiana conidia. We obtained and characterized the transcriptome of O. furnacalis larvae with special emphasis on immunity-related genes. 62,382 unigenes were assembled and 35,700 were annotated to known databases. Additionally, we performed quantitative reverse transcript (qRT)-PCR analysis to compare the gene expression profiles of B. bassiana-infected and non-infected O. furnacalis larvae. All these results give us an overview of gene expression profiles of O. furnacalis larvae response to B. bassiana, and provide a shortcut for identifying new immunity genes and useful information for studying the molecular basis of host-entomopathogenic fungus interaction.

Insect Rearing
Asian corn borer (O. furnacalis (Guenée)) was kindly gifted by Dr. Kanglai He from the Institute of Plant Protection, Chinese Academy of Agricultural Sciences. The larvae were reared on an artificial diet at 28uC under a relative humidity of 70-90% and a photoperiod of 16 h light and 8 h darkness [39].

B. bassiana Culture and Conidia Suspension Preparation
B. bassiana strain 252 was cultured on potato dextrose agar (PDA) plates at 25uC and 80% humidity. Conidia (spores) used for infection were harvested from 3-4 weeks old cultures by scraping the surface of the mycelia with sterile cell scrapers into sterile deionized water containing 0.1% Tween-80. Conidia were separated from other mycelial structures over a sterile funnel packed with autoclaved glass wool, washed twice with ddH 2 O by centrifugation at 4,000 rpm, counted and diluted to 2610 5 conidia/ml. Freshly prepared conidia were used for all experiments.

Immunization and Total RNA Extraction
Three microliter of diluted conidial suspension (2610 5 conidia/ ml) were injected into the haemocoel of O. furnacalis fifth instar day 0 larvae from the same batch. Injection of sterile deionized water was used as a control. After 10 h, each five larvae from challenged or control group were collected, and total RNA samples from the whole body were individually prepared using TRizol Reagent (TIANGEN, Beijing, China) following the manufacturer's instructions. Total RNA was dissolved in H 2 O, and RNA quantity was determined on a Nanodrop ND-2000 spectrophotometer (Nano-Drop products, Wilmington, DE, USA). RNA integrity was checked on Agilent 2100 BioAnalyzer (Agilent Technologies, Englewood, CO, USA).

Library Construction and Illumina Sequencing
Ten mg of total RNA equally from 5 larvae in each group was used to isolate mRNA using oligo(dT) magnetic beads. The cDNA library of each sample was constructed using NEBNextH mRNA Library Prep Reagent Set (NEB, Ipswich, MA, USA) following the manufacturer's protocols. Briefly, enriched poly(A) RNA of each sample was fragmented into 200-700 nt pieces with RNA Fragmentation Reagents. The cleaved RNA fragments were transcribed into the first-strand cDNA using random hexamerprimers, followed by second-strand cDNA synthesis. The resulting double-stranded cDNA (dsDNA) was purified with QiaQuick PCR extraction kit (Qiagen, Hilden, Germany) and resolved in EB buffer. The purified dsDNA was treated with T4 DNA Polymerase and T4 Polynucleotide Kinase for end-repairing and dA-tailing. After that, they were ligated to sequencing adaptors with barcode using T4 DNA ligase. Finally, fragments with around 200bplength were purified with QiaQuick GelPurify Kit (Qiagen, Hilden, Germany), and used as templates for PCR amplification to create the cDNA library. The library was paired-end sequenced using PE90 strategy on Illumina HiSeq TM 2000 (Illumina, San Diego, CA, USA) in the Beijing Genome Institute (Shenzhen, China). The challenged and control libraries were sequenced in one lane then raw-reads were sorted out by barcodes.

Assembly and Annotation of Transcriptomes
Raw reads from each library were filtered to remove low quality reads and the sequence reads containing adapters and poly-A/T tails. The resulting clean reads were assembled to produce unigenes using the short reads assembling program -Trinity [40]. The used parameters were as follow: min_glue = 2, V = 10, edge-thr = 0.05, min_kmer_cov = 2, path_reinforcement_distance = 80, and group_pairs_distance = 250. The other parameters were set as the default. Only those contigs with the length no shorter than 48 bp were used for the further assembly. The unigenes from the two samples were pooled together and clustered by TGI Clustering Tool [41]. The following parameters were used to ensure a high quality of assembly: a minimum of 95% identity, a minimum of 35 overlapping bases, a minimum of 35 scores and a maximum of 25 unmatched overhanging bases at sequence ends. The consensus cluster sequences and singletons make up the final unigene dataset.
For functional annotations, we firstly searched all unigene sequences against various protein databases such as Nr, Swiss-Prot, COG, and KEGG using BLASTX, and then searched nucleotide database Nt using BLASTN, with an E-value cut-off of 10 25 [42]. The BLAST results were used to extract coding region sequences (CDS) from the unigene sequences. The predicted CDS were further translated into peptide sequences. When a unigene happened to have no BLAST hits, ESTScan software [43] was used to determine the sequence direction. In addition, we performed the Gene Ontology (GO) annotations for each unigene using Blast2GO program, according to the GO associations from the BLASTX against the Nr database [44,45].

Identification of Differentially Expressed Genes
A mapping based on expression profile comparison was performed to examine transcript level change between control and treated sample. Clean reads of each library were mapped back to the common non-redundant dataset, respectively, using Bowtie allowing up to 3 base mismatches and a minimum length of 40 bp [46]. The expression of unigene was calculated with FPKM (Fragments per kb per million fragments) method [47]. P-values were calculated according to the hypergeometric test described by Audic and Claverie [48]. The calculated P-values were used to identify differentially expressed genes. The false discovery rate (FDR) was calculated according to the Benjamini and Hochberg algorithm [49], and further used to determine the P-value threshold in multiple testing. Genes with expression changes no less than two folds and FDRs less than 0.001 were considered as significant differentially expressed genes.
Identification and Sequence Analysis of Immunity-related Genes from O. furnacalis Transcriptome The available immunity-related gene sequences from other insect species were used as references to screen O. furnacalis transcriptome database obtained above. The referred insect species mainly included A. gambiae, A. mellifera, B. mori, D. melanogaster, M. sexta, T. castaneum, and T. molitor. The potential candidates of O. furnacalis immunity-related genes were confirmed by searching the BLASTX algorithm against the non-redundant NCBI nucleotide database using a cut-off E-value of 10 25 .
For the sequence analysis of putative immunity-related genes, the deduced protein domains were determined by using Pfam (http:// www.sanger.ac.uk/Software/Pfam/) and SMART (http://smart. embl.de/). Analysis of deduced amino acid sequences, including prediction of signal peptide, molecular weight and isoelectric point, was carried out in the EXPASY (Expert Protein Analysis System) proteomics server (http://www.expasy.org). Sequence comparisons and phylogenetic analysis were performed by MEGA5 software [50]. Phylogenetic trees were constructed by the neighbor-joining method, with statistical analysis by the bootstrap method, using 1000 repetitions. The phylogenetic trees were visualized in MEGA5 and colored in Adobe Illustrator (Adobe Systems).

Quantitative Reverse Transcript (qRT)-PCR Analysis
To validate expression profile from comparative transcriptomic analysis, we designed specific primers (Table S1) to perform qRT-PCR analysis for 24 immunity-related genes. The annealing temperatures of the primers were controlled at around 62uC. Total RNA from whole body was extracted independently from 3 biological replicates as described above. DNase I-treated RNA (1 mg) was converted into first-strand cDNA using TIANScript RT Kit (TIANGEN, Beijing, China). The cDNA products were diluted 5 fold for use as template. The qRT-PCR was performed on a Applied BiosystemsH 7500 Real-Time PCR System (Life Technologies, Grand Island, NY, USA) using the GoTaqH qPCR Master (Promega, Madison, WI, USA), according to the manufacturer's instructions. O. furnacalis ribosomal protein L8 (rpL8) was used as an internal standard to normalize the expression levels. The thermal cycling conditions for qRT-PCR were 95uC for 2 min followed by 40 cycles of 95uC for 15 s, 55uC for 30 s and 68uC for 30 s ending with a melting curve generation (60uC to 95uC in increment of 0.5uC very 5 s). The relative expression of genes was calculated using 2 2DDCt method [51].

Sequencing and Unigene Assembly
Parasitization of corn borer by B. bassiana provides a good system to study the interactions between insect hosts and entomopathogenic fungi, but the lack of genomic data of corn borer retards the relative progresses. In order to obtain detailed information about corn borer transcriptome, we subjected cDNA from O. furnacalis larvae with or without the infection of B. bassiana to Hiseq 2000 sequencing. A total of 57,411,104 and 57,669,432 raw reads were generated from water-injected (control) and B. bassiana-injected (treated) O. furnacalis libraries, respectively (Table  S2). After removal of adaptor sequences, ambiguous reads, and low-quality reads (Q20,20), control and treated libraries yielded 51,594,958 (SRA accession number SRX378863) and 52,437,534 (SRA accession number SRX378865) high-quality clean reads comprised of 4,643,546,220 nucleotides (4.64 Gb) and 4,719,378,060 nucleotides (4.72 Gb), respectively (Table S2). All high-quality reads were assembled de novo into 95,070 (control) and 96,561 (treated) contigs with a mean length of 371 and 352 nt, respectively (Table 1). Using paired-end reads and gap-filling, these contigs were further assembled into 66,004 (mean length 588 nt) and 71,723 (mean length 511 nt) unigenes (Table 1). These two unigenes-sets were pooled for further clustering, and finally revealed a common non-redundant dataset containing 62,382 unigenes with a mean length of 729 nt, which consist of 22,889 distinct clusters and 39,493 distinct singletons ( Table 1). The assembled sequences have been deposited in the NCBI Transcriptome Shotgun Assembly (TSA) Database with the title as BioProject: 228958TSA. The detailed information of each unigene, including gene ID, length, expression and functional annotation, was integrated in Table S3.

Gene Identification, Functional Annotation and Classification
All 62,382 unigenes were annotated by searching against the databases of Nr, Nt, Swiss-Prot, KEGG, COG, and GO. As shown in Table S4 (26,682,42.8%) was not annotated to the existing databases. It suggested that they were the potential sources of novel genes.
Functional classifications of all unigenes were determined by assigned to Gene Ontology (GO). GO is an international standardized gene functional classification system and covers three categories: biological process, cellular component, and molecular function. In our study, a total of 13,451 unigenes were assigned to one or more GO terms. Among these, 10,146 (75.4%) unigenes were grouped under the category of biological process, 7,510 (55.8%) under the category of cellular component, and 10,973 (81.6%) under the category of molecular function ( Figure  S2). The classification of GO terms was further performed at level 2 in each category. As occurred in the transcriptomes of other immune-activated insect larvae [36], the most abundant GO biological process categories were ''cellular processes'' (16.0%) and ''metabolic processes'' (12.6%). In this category, 437 unigenes comprise the subcategory of immune system process. In the cellular component category, ''cell'' (22.3%) and ''cell part'' (22.3%) represented the most abundant subcategories followed by ''organelle'' (14.9%). In the molecular function category, ''binding'' (41.4%) and ''catalytic activities'' (39.5%) were the most abundant ( Figure S2).

Identification of Differentially Expressed Genes in Response to B. bassiana Infection
To gain insight into the global transcriptional changes taking place in O. furnacalis larvae infected by B. bassiana conidia, we performed pairwise comparisons between water-injected and B.
bassiana conidia-injected libraries to identify the differentially expressed genes. The complete lists of differently expressed genes, including their GO and Nr annotations, are shown in Table S5. The results revealed that 13,890 unigenes exhibited significant changes after B. bassiana infection, including 5,843 up-regulated and 8,047 down-regulated unigenes. The detected fold changes (log2ratio) of gene expression ranged from 217.66 to 19.77.
The GO enrichment was performed to analyze the functions of all identified differently expressed genes. Among the 5,843 upregulated unigenes, 788, 582, and 910 ones were assigned to the categories of biological process, cellular component, and molecular function, respectively. Similarly, 8,047 down-regulated unigenes were also assigned to these three categories: 1,359 in ''biological process'', 943 in ''cellular component'', and 1561 in ''molecular function'' ( Figure 1 and Table S5). For both up-and downregulated unigenes, the top 2 most abundant subcategories in each GO category were as follows: ''cellular process'' and ''metabolic process'' in ''biological process'' cluster, ''cell'' and ''cell part'' in the ''cellular component'' cluster, and ''binding'' and ''catalytic activity'' in the ''molecular function'' cluster ( Figure 1). In the category of biological process, 42 up-regulated and 58 downregulated genes belonged to the term of immune system process.
To validate the data about differently expressed gene, qRT-PCR analysis was performed using specific primers for 24 immunity-related unigenes including 12 up-regulated, 4 downregulated and 8 unchanged unigenes. Data were presented as fold changes normalized to the rpL8 gene in the B. bassiana conidiainjected sample relative to the water-injected control sample ( Figure 2). Most tested unigenes exhibited the consistent expression trend in qRT-PCR analysis and the original differently expressed genes analysis ( Figure 2 and Table S5). It indicates that the results of gene expression profiling from transcriptome analysis are reliable. It is worth noting that the expression profiles of 4 genes including CL1725.Contig1 (OfCTL7), Unigene9842 (OfSpz-1A), unigene13709 (OfToll1), and CL997.Contig1 (OfPPO2) were different between qRT-PCR and transcriptome analysis. In the transcriptome analysis, Unigene9842 (OfSpz-1A) and CL997.Con-tig1 (OfPPO2) were slightly decreased while CL1725.Contig1 (OfCTL7) and unigene13709 (OfToll1) remained unchanged (Table  S5). In the qRT-PCR analysis, all four genes displayed a significant increase in transcript abundance after B. bassiana challenge ( Figure 2). The observed differences in gene expression might be caused by the difference in the accuracy of these two assay methods. The result from the qRT-PCR analysis could be more accurate because the differently expressed genes identified in this study were from the transcriptome assembly and mapping, but not from the special DGE analysis.

Identification of Immunity-related Genes
Insect innate immune response acts as a crucial part in defending against the infection of pathogens and parasites. In order to obtain the comprehensive perspective on the molecular biology of immune systems in O. furnacalis, we combined GO annotation and BLAST searches to identify the genes related with the cellular and humoral immune response. Based on molecular functions, we have divided the immune-related genes into 4 main groups: genes for signal recognition; genes involved in signal modulation and amplification; genes for signal transduction; and effector genes [52]. In total, we have identified 190 unigenes with high similarity to immunity-related genes, including 45 ones for signal recognition, 33 for signal modulation and amplification, 46 for signal transduction, and 66 for immune effectors (more details seen below). The deduced amino acid sequences of these 190 putative innate immunity genes are listed in Figure S4.
Genes for signal recognition. In insect innate immune response, recognition of nonself is the initial process. It is notably known that the recognition step is mediated by a group of proteins, known as pattern recognition proteins (PRPs), such as peptidoglycan recognition proteins (PGRPs), b-1,3-glucan recognition protein (bGRPs)/gram-negative binding proteins (GNBPs), Ctype lectins (CTLs), scavenger receptors (SCRs) and so on [24]. In O. furnacalis transcriptom, we totally identified at least 45 PRP transcripts including 10 PGRPs, 4 bGRPs, 14 CTLs, 9 SCRs, 2 hemocytins, 1 hemolin, 2 galectins, 1 dscam, 1 draper and 1 eater ( Table 2). PGRPs play a central role for the recognition of invading microorganisms in insect immunity by specifically binding to and hydrolyzing bacterial peptidoglycan [53,54]. The first PGRP was isolated from hemolymph of the silkworm, as a pattern recognition receptor to trigger prophenoloxidase (PPO) activation cascade [53]. All PGRP family members share at least one conserved PGRP domain, with similarity to bacteriophage T7 lysozyme, a zinc-dependent N-acetylmuramoyl-L-alanine amidase [55]. The most highly diversified PGRP homologues have been identified in Drosophila. Drosophila has 13 PGRP genes encoding 19 proteins, which are classified into short (S) and long (L) forms [55,56]. Among 19 Drosophila PGRPs, six (DmPGRP-SB1/2SB2/2 SC1a/2SC1b/2SC2/2LB) have amidase activity, and five (DmPGRP-SA/2SD/2LC/2LE/2LF) lack amidase activity but function as receptors to activate immune signaling pathways [55]. In this study, we identified 10 putative PGRP sequences and designated them as OfPGRP1-10. With the exception of OfPGRP2, the other 9 PGRP transcripts were predicted to be full-lengthed ( Table 2). Alignment of 10 putative O. furnacalis PGRPs with Drosophila PGRPs and T7 lysozyme indicated that the deduced amino acid sequences of OfPGRP1-3 and OfPGRP8-10 lack at least one of five active site residues essential for amidase activity in T7 lysozyme (H17, Y46, H122, K128 and C130, K128 is replaced by T in Drosophila PGRPs). It suggests that these six O. furnacalis PGRPs potentially act as receptors for peptidoglycan to initiate a signaling pathway while the left 4 ones (OfPGRP4-7) theoretically have amidase activity and might serve as an intracellular peptidoglycan scavenger. Bootstrap analysis reveals that OfPGRP1 is an ortholog of M. sexta PGRP-1 and B. mori PGRP-S1 which both have been verified to function as recognition receptors in the PPO activation cascade [53,54] (Figure 3). It suggests that OfPGRP1 may also act as a peptidoglycan receptor in activating PPO cascade upon the challenge of B. bassiana.
Moreover, the analysis of digital expression profile showed that 7 out of 10 identified PGRPs (OfPGRP1, 4-7, 9 and 10) were obviously up-regulated after the infection of B. bassiana, whereas the other 3 OfPGRPs remained unchanged ( Table 2 and Table S5). We randomly selected 4 PGRP genes (OfPGRP2, 6, 7 and 10) to analyze their transcript changes after the injection of B. bassiana using qRT-PCR methods. OfPGRP-6, -7 and -10 showed very high expression levels in B. bassiana-injected corn borer larvae, while OfPGRP2 mRNA level was consistent ( Figure 2). We speculated that O. furnacalis PGRPs might play different roles, and work in concert with each other to defend against the invasion of B. bassiana.
bGRP/GNBP belongs to another pattern recognition protein family. This family contains two functionally different proteins, one of which has a strong affinity to b-1,3-glucans of fungal cell walls (bGRP) and the other is dubbed as Gram-negative binding protein (GNBP) but binds to Gram-negative bacteria or Grampositive bacteria [57]. Since the first bGRP was identified in the PPO-activation system of B. mori [58], bGRPs have been identified in insects including Drosophila (3 genes) [57], Anopheles (7 genes) [27], Apis (2 genes) [29], Manduca (2 genes) [59,60], and Tribolium (3 genes) [30]. They all consist of a conserved N-terminal b-1,3glucan-recognition domain for the detection of pathogens or parasites, and a C-terminal glucanase-like domain with undefined function [61,62]. Three-dimensional structures of BmbGRP1 and DmGNBP3 further revealed that the N-terminal b-1,3-glucanrecognition domain actually adopts a b-sandwich structure formed by eight b-strands [62,63]. In this study, we identified 4 bGRP/ GNBP genes with predicted full length, and designated them as OfbGRP1-4. All OfbGRPs are assumed to be secreted proteins because they have putative signal peptides (Table 2). A comparison of the deduced amino acid sequences with Drosophila GNBP1-3 and Bombyx bGRP1 showed that OfbGRP1-3 contained the      OfbGRP1-3 also comprise eight conserved b-strands ( Figure 4A), suggesting their ability to bind to b-1,3-glucan. We performed the phylogenetic analysis for all OfbGRP1-4 and 42 bGRP sequences from other insect species to investigate the relationship between O. furnacalis bGRPs and others. As shown in Figure 4B, all 46 bGRPs were clustered into two groups, one containing the active catalytic residues for b-1,3-glucanase activity and the other containing no such residues. OfbGRP1-4 is presented as 1:1 orthologs to B. mori bGRP1-4. The orthologs for OfbGRP1, OfbGRP2 and OfbGRP3 are present in lepidopterans but not in dipterans ( Figure 4B). OfbGRP1 and OfbGRP3 are predicted to be paralogs and have expanded from a common ancestral gene after Lepidoptera had diverged ( Figure 4B). The digital expression profiles showed OfbGRP2 and OfbGRP4 mRNA level slightly increased while OfbGRP1 and OfbGRP3 were consistent in response to the infection of B. bassiana ( Table 2). The qRT-PCR assay showed OfbGRP1 mRNA level decreased and OfbGRP4 expression kept unchanged after challenge ( Figure 2). C-type lectin (CTL) is probably the largest lectin family. They are a group of soluble and membrane-bound proteins that associate with carbohydrates in a Ca 2+ -dependent manner [64]. Invertebrate CTLs participate in immune responses including PPO activation [65], hemocyte-mediated nodule formation and encapsulation [66], opsonization and microbial clearance [67]. CTLs have a characteristic carbohydrate-binding domain (CRD) with a well-defined structure stabilized by two or three pairs of disulfide bonds [65]. Insect CTLs normally consist of tandem CRDs. In this study, we identified 14 CTL genes, and designated them as OfCTL1-14 ( Table 2). All deduced O. furnacalis CTLs contain two consecutive CRDs, a characteristic for lepidopteran CTLs ( Figure 5). Although N-and C-terminal CRDs from Hyphantria cunea lectin are reported to have different sugar-binding specificities, the detailed function of these CRDs is not yet clarified [68]. We postulated that the two CRDs of O. furnacalis CTLs also have different sugar-binding specificity and bind to different microorganisms. Except for CRDs, other structural modules have been identified in some CTLs. For example, Drosophila Furrowed (a Drosopohila CTL) and Bombyx CTL2 have Sushi motifs and transmembrane domain in addition to CRD [31,69]. We only found CRDs in OfCTL1-14 sequences. There might be some unknown CTLs we did not identify in the transcriptome. OfCTL6 has 56% amino acid sequence identity to M. sexta IML-2 which binds to lipopolysaccharide to stimulate hemocyte encapsulation and melanization [64,65]. This suggests that OfCTL6 might also play an important role as recognition receptor during the early phase of microbial infection. Phylogenetic analyses indicated that OfCTL1-9 only clustered with lepidopteran CTLs, and OfCTL10-14 grouped with both lepidopteran and dipteran CTLs ( Figure 5). qRT-PCR assay revealed that the transcript level of OfCTL7 and OfCTL12 increased, OfCTL3 and OfCTL6 decreased, whereas OfCTL1 remained unchanged ( Figure 2).
In addition, we also identified other PRP genes from O. furnacalis transcriptome, including 9 scavenger receptors (SCRs), 2 hemocytins, 1 hemolin, 2 galectins, 1 dscam, 1 Draper and 1 Eater ( Table 2). The identification of a large number of PRP genes laid a good foundation for the further cloning and functional studies of PRPs in O. furnacalis.
Genes for extracellular signal modulation and amplification. After infectious parasites or pathogens are recognized by insects, the invasive signals are modulated by extracellular cascades. Similar to the coagulation pathway and complement system in human, insect plasma factors such as serine proteases (SPs) and serine protease homologs (SPHs) play critical roles in signal relaying/tuning and execution mechanisms [31,70]. SPs typically contains a catalytic triad consisting of His, Asp, and Ser residues, which are embedded in highly conserved sequence motif of TAAHC, DIAL and GDSGG, respectively [71]. SPHs lack proteolytic activity due to the substitution of the catalytic triad residues, but they can enhance specificities and catalytic activities of SPs [72][73][74]. SPs and SPHs constitute one of the largest protein families in insects [27,75,76]. We have identified 121 potential SP and SPH transcripts in O. furnacalis transcriptome. Given that the catalytic triad is vital to define a SP or SPH, we manually checked the deduced amino acid sequences of all 121 transcripts, and only kept those containing all three sequence motifs. Finally we obtained 47 SP (designated as OfSP1-OfSP47) and 14 SPH (designated as OfSPH1-OfSPH14) transcripts. Fifty-six genes are predicted to have full length, 18 of which encode polypeptides with a SP or SP-like domain and other structural modules. These include 13 SPs (SP1-SP5, SP7, SP8, SP10, SP12-SP14, SP17 and SP37) and 3 SPHs (SPH8-SPH10) which contain one or more regulatory clip domains (Table 2), one SP (SP40) which contains a CUB domain, one SP (SP46) which contains five low density lipoprotein receptor A repeats (LDLa) domains and two complement control protein (CCP) domains. The clip domain is an important structural unit in which six conserved cysteine residues form three disulfide bonds [25,77]. In arthropods, clip-domain SPs, and occasionally clip-domain SPHs, are involved in many immune signaling pathways, such as melanization cascade and Toll pathway [31,71,76]. In 13 SPs with clip domain, SP1 and SP13 mediated the immune responses of corn borer against B. bassiana by participating in the PPO activation cascade (submitted to Amino Acids). The other 11 clip-SPs are still under the investigation. Moreover, it is noteworthy that O. furnacalis SP46, with a large size and complex domain structure, is most similar to M. sexta HP14 [78] and T. molitor MSP [21] which both function as an initial enzyme to be recruited into the recognition complex in PPO activation cascade. We, therefore, inferred that SP46 also acts as the first enzyme in a serine protease pathway. Phylogenetic analysis indicated that O. furnacalis clip-domain proteins are divided into four subfamilies ( Figure S5). The subfamily A is composed of SPHs solely while subfamilies B, C and D comprise SPs mainly. The four groups of SP-related genes may represent lineages derived from ancient evolutionary events since similar subfamilies also existed in Anopheles [27], Drosophila [75], and Tribolium [30].
Genes involved in signal transduction. After the invasive signals from the infectious microorganisms are recognized and modulated, the signal transduction pathways will be initiated to produce the effector molecules. Four signal transduction pathways, Toll, Imd, JNK, and JAK/STAT are known to be involved in insect immunity [29]. Toll and Imd pathways are more important in sensing microbes. The Toll pathway is primarily involved in the defense against fungi and Gram-positive bacteria with lysine-type peptidoglycans (Lys-type PGNs) in their cell walls, while the Imd pathway responds to Gram-negative bacteria and some Grampositive bacteria with meso-di-aminopimelic acid-type peptidoglycans (DAP-type PGNs) [1]. In this study, we have identified at least 46 transcripts for signal transduction from O. furnacalis transcriptome, which nearly included all known components in signal transduction pathways (Table 2).
Spä tzle is a ligand for the Toll receptor and activates the Toll signaling pathway [11]. Spätzle is present in Anopheles (6 genes) [27], Drosophila (6 genes) [26], and Bombyx (3 genes) [31]. Drosophila Spz1, Bombyx Spz1, and Manduca Spz1 all have been demonstrated to participate in immunity [89][90][91]. ProSpä tzle is consisted of an unstructured pro-domain and a C-terminal fragment that adopts a cysteine-knot structure similar to that of mammalian neurotrophins [92]. Six Spätzle genes have been identified from the transcriptome, with the tentative name as OfSpz1A, 1B, and 3-6. Only OfSpz1A, 1B, and -5 contains the complete open reading frame (Table 2). To assess the relationship between O. furnacalis Spä tzle and other insect Spä tzle proteins, we performed a phylogenetic analysis by aligning the known homologous cysteine-knot domain sequences from different insect species. The phylogenetic tree ( Figure 6) suggests that all Spä tzle homologs can be assigned to a 1:1 orthologous group with one of the Drosophila Spä tzle gene products (Spz1-Spz6). No Spä tzle-2 ortholog was identified in O. furnacalis transcriptome. A possible reason is that Spä tzle-2 gene is missing in O. furnacalis because of the evolutionary event. The other reason with higher possibility is that the transcript level of Spä tzle-2 is low and it is not captured in the RNA-seq. As shown in Figure 6, the bootstrap value at the clade including Spz-1 is lower than in the other clades containing Spz-2 to Spz-6, indicating a lower degree of sequence conservation in Spä tzle-1. It is worth noting that there are two possible O. furnacalis Spä tzle-1 variants (OfSpz-1A and OfSpz-1B) in Spä tzle-1 branch, which only share 33% identity in amino acid sequences. OfSpz-1A is more similar to other insect Spä tzle-1, which had 56% sequence similarities with Manduca Spä tzle-1. qRT-PCR analysis indicated that the expression of both OfSpz-1A and OfSpz-1B were induced by the B. bassiana conidia. It suggests that certain genes in Toll pathway are involved in innate immunity response against B. bassiana.
Toll receptor plays a crucial role in insect innate immune response and acts as the signal transducer in the Toll pathway [93]. It has been identified in many insect species, including the orthoptera, hymenoptera, coleoptera, lepidoptera and diptera. Toll-like receptors (TLR) have also been detected in mammals, suggesting that Toll and TLR are evolutionarily conserved throughout the insects and mammals [94]. Insect Toll receptors and mammalian TLRs are all type I membrane proteins with an ectodomain consisting of leucine-rich repeats (LRRs), a transmembrane domain, and an intracellular Toll-interleukin homolog domain (TIR) that can transduce signals [95]. In this study, we identified 12 genes coding Toll receptors in O. furnacalis transcriptome datasets. These genes were designated as OfToll1-12. Only OfToll1-4 and OfToll8 genes contain the full-length encoding sequences (Table 2). Domain prediction with SMART revealed that among these 5 full-lengthed Toll-like genes, OfToll-1 through -3 is composed of the extracellular LRR, transmembrane and cytoplasmic TIR domains, while OfToll4 and OfToll8 only contains the extracellular LRRs and transmembrane domains (Figure 7). Among five incomplete O. furnacalis Toll-like genes with stop codon at 39-end (OfToll5-7, -9, and -11), cytoplasmic TIR domain is only found in OfToll-5, OfToll-7, and OfToll-11 ( Figure 7). We cannot conclude whether TIR domain is present in OfToll-10 and OfToll-12 because both 59-and 39-end of the cDNA sequences of these two genes are incomplete (Table 2 and Figure 7). On the other hand, the TIR domain has a more reliable determination of phylogeny than the extracellular LRR region [96]. Therefore, we were unable to conduct the convincible phylogenetic analysis for identified O. furnacalis Tolls. The efforts are being made to obtain the complete coding sequences for all potential O. furnacalis Toll genes. We investigated Toll gene expression after B. bassiana infection. Entopathogenetic fungus B. bassiana significantly increased the mRNA level of Toll-1, but decreased the transcript level of Toll-5 ( Figure 2).
In Toll signaling pathway, the cysteine-knot part of Spä tzle binds to the ectodomain of the Toll receptor and thereby triggers an intracellular signaling cascade including dMyD88, Tollip, Tube, Pelle, Pellino, TRAF2, and finally results in the degradation of inhibitor protein Cactus and the nuclear import of a rel family transcription factor Dorsal/Dif [9,97]. Contrary to the ligandreceptor diversification, components of the intracellular cascade appear to be highly conserved in insects [30]. In O. furnacalis transcriptome dataset, we have identified dMyD88, Tollip, Tube, Pelle, Pellino, TRAF2 with 1:1 orthologs (Table 2). We postulated that a similar protein complex also forms in O. furnacalis, which will lead to the release of a Rel transcription factor (CL7451.Contig1) from a cactus-like protein (CL1201.Contig1), allowing the Rel molecular to translocate into the nucleus and activate the expression of effector genes such as antimicrobial peptides ( Figure  S7).
Imd pathway is also composed of multiple molecules for the signal transduction. The involved intercellular components include Imd, Fas-associated death domain protein (FADD), Dredd, IAP2, transforming growth factor b activated kinase (TAK1), Tab2, Ubc13, an inhibitor of nuclease factor kB kinase subunits b and c (IKKb and IKKc), and a transcription factor Relish [98,99]. Orthologs of all these components were identified from O. furnacalis transcriptome ( Table 2 and Figure S7). It strongly suggests that Imd-mediated immunity is conserved in Asian corn borer. Further experiments are required to verify the suggested role of each component in the immune response of O. furnacalis.
Genes for immune responsive effectors. Effector genes are induced in some specific tissues, such as fat bodies and hemocytes, following successive immune processes of signal recognition, modulation, and transduction. Part of the synthesized effector molecules are released into the hemolymph to play direct roles in phenoloxidase-dependent melanization, elimination of infectious microorganisms, apoptosis and other immune-related mechanism.
Prophenoloxidases (PPOs) are copper-containing enzymes. They are synthesized as inactive zymogens and activated by cleavage after residue arginine at approximately residue 50 [100,101]. A serine protease cascade contributes to the PPO activation. Active phenoloxidase (PO) catalyzed the hydroxylation of monophenols to o-diphenols and the oxidation of o-diphenols to quinones. Quinones are involved in microbial killing, melanin synthesis, sequestration of parasites or pathogens, and wound healing [7]. We have identified three PPO transcripts, and designated them as OfPPO1-3 ( Table 2). The total number of PPO genes in different insect species did not change significantly, except for the mosquito Anopheles: 3, 9, 1, 2, 2, 2 PPO genes have been identified from Drosophila, Anopheles, Apis, Bombyx, Manduca, and Tribolium, respectively [31]. Among three identified PPO fragments, only OfPPO2 sequence was confirmed in other studies [102]. OfPPO2 is also the only one with complete coding region. Signal sequences of OfPPO1-3 are not confirmed, consistent with the situation in other insects. It suggests that PPOs are released from cells into the hemolymph by cell rupture [100]. In addition, PPO subunits form heterodimer in Bombyx (BomPPO1 and 2) and Manduca (MsPPO1 and 2), but form homodimer in Drosophila (DmDox-A3) [103]. We infer that Ostrinia PPO subunits can also form such heterodimer or homodimer. Phylogenetic analysis showed that O. furnacalis PPOs cluster well with other lepidopteran PPOs. The 1:1:1 orthologous relationship existed between OfPPO1, BmPPO1, and MsPPO1, as well as between OfPPO2, BmPPO2, and MsPPO2. It suggested that an ancestral PPO gene was duplicated after lepidoptera had diverged (Figure 8). qRT-  groups: cecropins, insect defensins, lysozymes, proline-rich proteins, and glycine-rich proteins such as attacin [104,105]. Analysis of the O. furnacalis transcriptome revealed a large number of unigenes with homology to the major families of AMPs such as lysozymes, cecropins, moricins, lebocins, gloverins, defensins, and attacins (Table 2).
Lysozymes catalyze the hydrolysis of the b-1,4-glycosidic bond between N-acetylglucosamine and N-acetylmuramic acids of Figure 9. Phylogenetic analysis of lysozymes. The names of lysozyme genes used in the analysis were shown as scientific name of species followed by GenBank accession number of this specific gene. The Ostrinia lysozymes are marked in red. The branches specific for invertebrate c-type, vertebrate c-type, i-type, and g-type lysozymes are shaded in yellow, blue, green, and orange, respectively. For explanation of the arrows see Fig. 6. doi:10.1371/journal.pone.0086436.g009 peptidoglycans in bacterial cell walls [106]. They are divided into five major types including the c-type (chicken type), g-type (goose type), i-type (invertebrates), phage type, bacteria type, and plant type [107]. In this study, we identified six putative lysozyme genes (OfLys1-Lys6), and two putative lysozyme-like genes (OfLys-L1 and OfLys-L2). OfLys-L1 and OfurLys-L2 have an amino acid replaced at one or both of catalytic amino acid position (glutamic acid and aspartic acid), and, therefore, may lose muramidase activity [108]. The phylogenetic analysis revealed that the c-, g-and i-type lysozymes form three independent clusters, respectively (Figure 9). Among c-type cluster, the invertebrate and vertebrate formed a subgroup respectively. OfLys1-OfLys4 is included in invertebrate c-type lysozyme subgroup. OfLys5 and OfLys6 are clustered with i-type lysozymes with a high bootstrap support. We randomly selected OfLys2 for qRT-PCR assay. The result indicated that the expression of OfLys2 was significantly up-regulated upon the challenge of B. Bassiana (Figure 2).
Cecropins represent another group of linear and amphipathic peptides which adopt a-helical structure. The first cecropin was discovered in Hyalophora cecropia [109]. They are widely present and active against Gram-negative and Gram-positive bacteria and fungi [110]. Both Drosophila and Anopheles contain 4 cecropin genes and Apis has none [27,29]. However, up to 13 cecropin genes was detected in Bombyx genome [31]. In this study, we have identified six putative cecropin genes, designated as OfCec1-OfCec6 (Table 2). All identified Ostrinia cecropin genes contain complete coding regions, and their deduced amino acid sequences have a predicted signal peptide with 22 residues. Alignment of the conceptive protein sequences of OfCec1-6 showed that they have high similarities. For example, OfCec1 and OfCec4 share 82% identical amino acid residues, OfCec2 and OfCec3 have 85% identity in amino acid sequences ( Figure 10A). Similar to the situation observed in the cecropin family in Spodoptera exigua [34], it is difficult to clearly distinguish that the sequence variations represent either different genes or different alleles (or alternative transcripts) from the same gene. Phylogenetic analysis showed that Drosophila, Anopheles and lepidopteran insects form cecropin gene clusters ( Figure 10B). qRT-PCR assay for OfCec3 showed that the expression level of this gene significantly increased in response to B. bassiana infection (Figure 2). Other than lysozyme and cecropin genes, we have identified several other unigenes encoding putative AMPs in O. furnacalis transcriptome dataset, including four moricins, three lebocins, two gloverins, one defensin and one attacin ( Table 2). Expression profile analysis with qRT-PCR methods showed that all tested antimicrobial peptide genes were significantly induced in response to the B. bassiana injection, with the exception of defensin remaining unchanged ( Figure 2).
In addition, oxygen-derived free radicals and gaseous radical nitric oxide are also potent immune effector molecules [111]. Therefore, genes encoding the enzymes for the conversion of reactive oxygen species (ROS) and reactive nitrogen species (RNS) are somewhat classified as immune responsive effector genes. These enzymes include peroxidase (Pox), glutathione oxidase (GTX), superoxide dismutase (SOD), catalase, thioredoxin, thioredoxin reductase, peroxiredoxin, nitric oxide synthase (NOS), NADPH oxidase (NOX) and so on [30]. We have identified some of these genes in O. furnacalis transcriptome, including 17 Pox, 6 SOD, 1 catalase, 4 thioredoxin, 1 thioredoxin reductase, 5 peroxiredoxin, 2 NOS, and 1 NOX (Table 2). This suggests that local production of free radicals might be also a critical component of the immune response in Asian corn borers.

Conclusions
In summary, we sequenced and characterized the transcriptome from water-injected and B. bassiana-injected Asian corn borers. The transcriptome datasets obtained in this study make a significant contribution to a comprehensive sequence source for future O. furnacalis study, especially under the situation where its genomic information is currently unavailable. The explored immunity-related genes constitute an integrated picture of the immune network, which provides the valuable clues for a better understanding of the immune processes in O. furnacalis against B. bassiana. Immune constituent genes involved in signal recognition, modulation, transduction, and effector mechanisms have been identified and analyzed from the transcriptome. These immune repertoire genes appear to be evolutionarily conserved to different extent, and have various transcriptal profiles in response to the infection of B. bassiana. Functional analyses are necessary to verify our predictions. Nevertheless, the framework of information presented in this study should help clarify immune functions in an important agriculture pest and further understand the complex interaction between the insect pest and its entomopathogenic fungus.