Parasitization by Cotesia chilonis Influences Gene Expression in Fatbody and Hemocytes of Chilo suppressalis

Background During oviposition many parasitoid wasps inject various factors, such as polydnaviruses (PDVs), along with eggs that manipulate the physiology and development of their hosts. These manipulations are thought to benefit the parasites. However, the detailed mechanisms of insect host-parasitoid interactions are not fully understood at the molecular level. Based on recent findings that some parasitoids influence gene expression in their hosts, we posed the hypothesis that parasitization by a braconid wasp, Cotesia chilonis, influences the expression of genes responsible for development, metabolism and immune functions in the fatbody and hemocytes of its host, Chilo suppressalis. Methodology/Principal Findings We obtained 39,344,452 reads, which were assembled into 146,770 scaffolds, and 76,016 unigenes. Parasitization impacted gene expression in fatbody and hemocytes. Of these, 8096 fatbody or 5743 hemocyte unigenes were down-regulated, and 2572 fatbody or 1452 hemocyte unigenes were up-regulated. Gene ontology data showed that the majority of the differentially expressed genes are involved in enzyme-regulated activity, binding, transcription regulator activity and catalytic activity. qPCR results show that most anti-microbial peptide transcription levels were up-regulated after parasitization. Expression of bracovirus genes was detected in parasitized larvae with 19 unique sequences identified from six PDV gene families including ankyrin, CrV1 protein, cystatin, early-expressed (EP) proteins, lectin, and protein tyrosine phosphatase. Conclusions The current study supports our hypothesis that parasitization influences the expression of fatbody and hemocyte genes in the host, C. suppressalis. The general view is that manipulation of host metabolism and immunity benefits the development and emergence of the parasitoid offsprings. The accepted beneficial mechanisms include the direct impact of parasitoid-associated virulence factors such as venom and polydnavirus on host tissues (such as cell damage) and, more deeply, the ability of these factors to influence gene expression. We infer that insect parasitoids generally manipulate their environments, the internal milieu of their hosts.


Introduction
Parasitoid wasps of the order Hymenoptera develop as parasites of other arthropods during their larval stages, giving rise to freeliving adults. They are valued biological control agents for various insect pests [1]. Endoparasitoid wasps (whose larvae develop inside, rather than, on their hosts) introduce substances into their hosts during oviposition, including venom, polydnaviruses (PDVs), ovary fluids, and other maternal factors; these materials act to ensure successful development of their progeny [2]. These factors influence host behavior [3], metabolism, development [4], endocrine system activity [5] and immune defense reactions [1,6,7]. There are over 100,000 host-parasitoid systems and most of them are shaped by differing selective forces [6]. These coevolved systems have produced an unknown, but large number of variations on the broad theme of molecular host-parasitoid interactions. Only a few of these relationships have been deeply investigated, and much more knowledge is required to generate broad principles of molecular parasitoid-host systems.
The rice stem borer, Chilo suppressalis (Walker) (Lepidoptera: Crambidae) is a destructive rice pest in China and other Asian countries. It is responsible for severe crop loss every year, especially in China because of changes in rice cultivation and the popularization of hybrid rice. Hybrid strains are more susceptible to insect damage than other rice releases [8]. C. suppressalis has developed resistance to many groups of chemical insecticides [9,10] and the estimated cost of controlling this pest is around 1 billion yuan annually [11]. Crop damage and high resistance emphasizes the urgency for developing innovative control measures and resistance management strategies. Parasitoids or parasitoid-produced regulatory molecules have the potential to improve conventional pest control strategies in ways that supports sustainable agriculture [12]. Cotesia chilonis (Matsumura) (Hymenoptera: Braconidae), mainly distributed in southeastern and eastern parts of Asia, is the major endoparasitoid of C. suppressalis larvae [13]. C. chilonis injects venom, PDV and teratocytes as major parasitoid-associated factors while ovipositing into hosts [14]. The injected virus is in the genus Bracovirus (BV) (Family: Polydnaviridae) similar to Cotesia vestalis [15]. The biological characteristics of C. chilonis and its effects on the immune response of C. suppressalis larvae has been preliminarily investigated [16]. When dissecting the parasitized hosts, we found that the egg matured in 2 d, and larvae seemed to have three instars, the first two instar ones molted inside the host, and the third instar ones emerged from the host to spin a cocoon. The first, second, and third instar lasted 2, 3, and 1 d at 2561uC and 60,65% relative humidity, respectively. The pupae develop for 3 d. After parasitization by C. chilonis, total amount of food consumption of host larvae, compared with non-parasitized larvae, reduced by 36.75%. During the parasitization, the host development rate was restrained and the times of host-mounting become less, and the host larvae could not develop into pupae stage [17]. Parasitization by C. chilonis may also result in some regular changes of immunity of its host C. suppressalis [18]. For example, total number of hemocytes in parasitized larvae became significantly higher than that of non-parasitized (n.p.) control from 1 day postparasitization (p.p.) [18].
Parasitoid wasps have evolved an array of mechanisms to regulate the host's physiology and biochemistry in a way that creates a microenvironment for successful development [6]. Previous studies have concentrated mainly on individual or small defined groups of host genes to explore their functions or differential expression following parasitization [19][20][21]. Only a few studies report on large-scale approaches to understanding the global impacts of parasitization on hosts at the genome level. Using suppression subtractive hybridization, Fang et al. [22] found that Pteromalus puparum venom treatments led to reductions in expression of a large number of immune-related genes in the lepidopteran host Pieris rapae. Gene expression changes in flour moth Ephestia kuehniella caterpillars after parasitization by the endoparasitic wasp Venturia canescens were analyzed using cDNAamplified fragment length polymorphisms, which demonstrated that expression of 13 transcripts in parasitized hosts were suppressed by the wasp [23]. Deep sequencing-based transcriptome analysis of Plutella xylostella larvae parasitized by Diadegma semiclausum also indicated that parasitization had significant impacts on expression levels of 928 identified insect host transcripts [12].
In the present study, we used the Illumina sequencing technology to explore the C. suppressalis gene expression changes induced by C. chilonis parasitization. We first obtained and characterized the transcriptome of C. suppressalis larvae parasitized by C. chilonis. A systematic bioinformatics strategy was engaged to functional annotation of the transcriptome. Additionally, we constructed four RNA-seq (quantification) and compared the accumulation of transcription products of fat body and hemocytes in non-parasitization (n.p.) versus post-parasitization (p.p.) hosts, C. suppressalis. The results give us a comprehensive view of global gene expression profiles of two immune-related tissues of host response to parasitization, and establish a sound foundation for future molecular studies based on high throughput sequence data.

Materials and Methods
Insect rearing, parasitization and RNA isolation C. suppressalis were reared on artificial diet [24]. The wasps, C. chilonis, were reared on host larvae. Both species were maintained at 2561uC under natural photoperiod and relative humidity approximately 80%. To obtain material for sequencing, 100 larvae with the age of day 2 (4 th instar) were exposed to a mated female wasp until parasitization was observed. Individual parasitized larvae were maintained on artificial diet under the conditions described until tissue samples were prepared.
Larvae of C. suppressalis were surface-sterilized with 70% ethanol. Hemocytes were prepared by puncturing a proleg and allowing hemolymph to freely drip into insect Grace's medium (1:10, v/v; Invitrogen, Carlsbad, CA) in 1.5 ml chilled Eppendorf tubes and centrifuged at 200 6 g for 10 min at 4uC After centrifugation, plasma was discarded and hemocytes were used for total RNA extraction. The fat bodies were removed from the remaining cadaver under a stereomicroscope and transferred into phosphate-buffered saline (NaCl 137 mM, KCl 2.7 mM, Na 2 HPO 4 10 mM, KH 2 PO 4 2 mM, pH 7.2 , 7.4) in 1.5 ml Eppendorf tubes. Total RNA samples were extracted using TRIZOL Reagent (Invitrogen) following the manufacturer's instructions and stored in -80uC. RNA sample concentrations were determined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). Integrity was ensured through analysis on a 1.5% (w/v) agarose gel.

Transcriptome analysis library preparation and sequencing
The previous of our work showed that the immune indices like hemocyte spreading rate, mortality, phagocytic rate, encapsulation index and phenoloxidase activity were all significantly changed after parasitism in 0.5 to 2 days [18]. Besides this, immature development of C. chilonis was studied by dissecting parasitized hosts in the laboratory at 2561uC and 60 -65% RH. When dissecting the parasitized hosts, we found that the egg matured in 2 d [14]. Hence, we selected 6, 12, 24 and 48 hrs p.p. based on the influences of C. chilonis development on host immunity [18]. The purpose of these four time intervals was to obtain a comprehensive sampling of transcripts, some of which would have been missed if tissues were collected at a single time point. Cs-FB and Cs-HC RNA was prepared at the same time as fat body and hemocytes from parasitized larvae (PCs-FB; PCs-HC) from day 2, 4th instar naïve larvae (100). To obtain complete gene expression information, a pooled RNA sample including sixteen RNA samples composed of four time points (6, 12, 24 and 48 h) of four treatments (PCs-FB, PCs-HC, Cs-FB and Cs-HC) was used for transcriptome analysis.
The cDNA library was prepared according to the Illumina manufacturer's instructions. Briefly, oligo (dT) beads were used to isolate poly(A) mRNA from total RNA (pooled RNA of control and experimental fat body and hemocytes). Short mRNA fragments were created by adding fragmentation buffer. Then, first and second strand cDNA were synthesized from cleaved RNA fragments. Short fragments were purified with QiaQuick PCR extraction kits (Qiagen, Hilden, Germany) and resolved with EB buffer for end reparation and adding poly(A). The short fragments were connected to sequencing adapters. Following agarose gel electrophoresis, suitable fragments were selected for PCR ampli-fication as templates. The library was sequenced using Illumina HiSeq TM 2000 (Illumina, Inc, San Diego, CA) at Beijing Genomics Institute (BGI)-Shenzhen, China (http://www. genomics.cn).

De novo transcriptome assembly and unigene annotation
The raw reads from the images and quality value calculation were performed by the Illumina data processing pipeline (version 1.6). After removal of low quality reads, clean reads were assembled into sequence contigs, scaffolds, and unigenes using the short reads assembling program SOAPdenovo [25]. All raw sequencing data have been deposited in NCBI's Short Read Archive (SRA) database (http://www.ncbi.nlm.nih.gov/sra) under accession number: SRR651040. The Transcriptome Shotgun Assembly (TSA) project has been deposited at DDBJ/EMBL/ GenBank under the accession number: GAJS00000000.
The unigenes were used for BLAST search and annotation against NR database and Swissprot database with an E-value cut of E-value 25 . Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology (KO) annotations of the unigenes were determined using Blast2GO and Blastall software [26].

RNA-seq (quantification) library preparation, sequencing and alignment with references
Total RNAs of four time points were mixed equally to create one library. Therefore, four RNA-seq libraries including, PCs-FB, PCs-HC, Cs-FB and Cs-HC, were prepared. The RNA-seq sequencing method was the same with transcriptome analysis (transcriptome analysis library preparation and sequencing). Briefly, after filter procedures, we obtained the clean reads, which were the basis of all following analysis. For the BGI bioinformatics pipeline, clean reads from each library were separately mapped against the reference set of assembled transcripts using SOAPaligner/soap2 [25]. Mismatches of no more than 2 bases were allowed in the alignment.

Gene expression level and differentially expressed genes identification
Gene expression levels were calculated using Reads Per Kilobase per Million (RPKM) mapped reads [27]. If there was more than one transcript for a gene, the longest one was used to calculate its expression level and coverage. Thus, the output for each dataset can be directly compared as the number of mapped reads per dataset and transcript size has been taken into account.
The correlation of the detected count numbers between parallel libraries were assessed statistically by calculating the Pearson correlation. False discovery rate (FDR) was used to determine differentially expressed genes [28]. Assume that we have picked out R differentially expressed genes in which S genes show differential expression and the other V genes are false positives. If the error ratio Q = V/R must remain below a cutoff (1%), FDR should not exceed 0.01. In this research, P # 0.01, FDR # 0.001 and the absolute value of log 2 Ratio $ 1 were used as threshold values to identify differentially expressed genes [29].

Quantitative real-time PCR (qRT-PCR) validation
Total RNA was extracted as described for RNA-seq library preparation and sequencing. Following DNAse I (RQ1 RNasefree DNase: Promega) treatment, total RNA (1mg) was used for cDNA synthesis with ReverTra Ace qPCR RT kits (Toyobo, Osaka, Japan).. Quantitative RT-PCR (qPCR) reactions (20 ml) were performed in triplicate using SsoFast EvaGreen Supermix with low ROX (BioRad) in a 7500 Real Time PCR System (Applied Biosystems by Life Technologies). The qPCR reaction consisted of 2 ml of diluted cDNA (10 ng) and 1 mM of each primer, which were selected for at least 90% amplification efficiency. The PCR reactions were programmed at 95uC for 30 sec; 40 cycles of 95uC for 5 sec, 60uC for 34 sec, followed by melting curve analysis for quality control (60uC to 95uC). No primer dimer was detected in the melting curves. The data were analyzed using the comparative Ct (ddCt) method [30], and the endogenous 18S rRNA reference gene [31] was used for normalization. At least three replicates were tested per sample.
We performed another experiment to record gene expression levels at 6, 12, 24 and 48 h p.p. for a selected group of genes. For each time point, three independent groups of 30 control larvae and three independent groups of 30 parasitized 4 th instar larvae were processed for RNA extraction.

Results and Discussion
Illumina sequencing and reads assembly Illumina sequencing resulted in 39,344,452 raw reads, corresponding to an accumulated length of 3,541,000,680 bp ( Table 1). The raw reads were assembled into 1,028,924 contigs with a mean length of 127 bp. Using paired end-joining and gap-filling, these contigs were further assembled into 146,770 scaffolds with a mean length of 275 bp. Scaffold sequences were assembled into clusters using TGI software. We obtained 76,016 unigenes with a mean length of 440 bp. The lengths of 18,462 unigenes were $ 500 bp and the lengths of the remaining 57,554 unigenes (75% of the total) were between 100 to 500 bp (Figure 1), similar to other insect transcriptome projects using this technology [32,33].

Annotation of predicted proteins, GO and COG classification
For functional annotation, the 76,016 unigenes were searched using BLASTx, with a threshold of E value , 10 25 , against four public databases (NCBI non-redundant (nr) database, the Swiss-Prot protein database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, and the Clusters of Orthologous Groups (COG) of proteins database. The E-value distribution of the top hits in the nr database showed that 24% of the mapped sequences have strong homology (less than 1.0E 249 ) and 76% of homolog sequences ranged between 1.0E 25 to 1.0E 249 ( Figure  2A). The similarity distribution has a comparable pattern with 21% of the sequences having similarity higher than 80%, while 79% of the hits have similarities ranging from 28% to 80% ( Figure  2B). The results are similar to transcriptome analyses of other insect species using this technology [34,35]. The species distribution of the best match result for each sequence showed that 40% of the C. suppressalis sequences match with sequences from the Drosophila species, while very low proportion (2%) of them have matches to Bombyx mori ( Figure 2C). One reason for the higher number of hits against the fruit fly genome is that approximately ten times more Drosophila genes than B. mori genes are deposited in databases.
In total, 11,886 unigenes were assigned GO terms based on Blast2GO [26] and WEGO [27] software. In each of the 3 main categories of GO classification, biological process (cell process dominates), cellular component (cell part dominates), and molecular function (binding dominates), show the analyzed tissues were most likely undergoing rapid growth and extensive metabolic activities. We did not find genes representing other clusters. We registered a high-percentage of genes from categories of ''metabolic process'', ''biological regulation'' and ''catalytic activity'' and only a few genes from terms ''synapse part'' and ''antioxidant activity'' ( Figure S1). We assigned 14,809 unigenes to COG clusters ( Figure S2). Among the 25 COG categories, the cluster for ''General function prediction'' represents the largest group (2587, 17.5%) followed by ''Replication, recombination and repair'' (1438, 9.7%) and 'Translation, ribosomal structure and biogenesis' (1219, 8.2%). The category of ''secondary metabolites biosynthesis, transport and catabolism'' (414, 2.8%) was particularly important because of the importance of secondary insecticide metabolites in insects. The most abundant sequences in this category are cytochrome P450 monooxygenases.

Statistics of RNA-seq (quantification) and differential gene expression
To characterize the gene expression profiles of fatbody and hemocytes in parasitized C. suppressalis by C. chilonis, four RNA-seq (quantification) libraries were constructed and sequenced. We generated 12,052,737 reads from control fat body (Cs-FB), 12,361,322 from parasitized fat body (PCs-FB), 12,466,924 from control hemocytes (Cs-HC) and 11,471,001 from parasitized hemocytes (PCs-HC) ( Table 2). These reads were mapped with reference sequences. Our data analyses indicate that parasitism has a significant impact on the gene expression profile of larval fatbody and hemocytes. For fatbody, 10,668 unigenes were differentially expressed after parasitization, with 2,572 (24%) upregulated and 8,096 (76%) down-regulated. For hemocytes, 7,195 unigenes were differentially expressed after parasitization, with 1,452 (20%) up-regulated and 5,743 (80%) down-regulated ( Figure  3). It can be shown that only 14% transcripts of C. suppressalis were differentially expressed after parasitization. It indicated that parasitization alter the abundance of a relatively low proportion of C. suppressalis transcripts in fat body and hemocytes

GO analysis of differentially expressed unigenes
Most of the differentially expressed transcripts (DETs) for the GO terms, molecular function and biological process, were downregulated except antioxidant activity ( Figure 4). This finding differs from the analysis of P. xylostella parasitized by D. semiclausum because most of the DETs were up-regulated [12]. One reason may be that this is a species-specific response and another reason may be that different PDV genera are associated with these two parasitoid wasps. Ichnovirus (IV) is associated with D. semiclausum and PDV with C. chilonis belongs to BV. Although viruses in these two genera have similar immunosuppressive and developmental effects on parasitized hosts, they differ morphologically and their encapsidated genomes largely encode different genes [15,36].

Transcripts related to immunity
Parasitism exerted significant impact on the transcriptome profile of fatbody and hemocytes. Among the changed unigenes, those related to immunity, development and metabolism are displayed in Table 3 and 4. These transcripts are most relevant to parasitism.
There are two pathways for pathogen recognition and signal transduction, a PRR-SP system in insect plasma (e.g., spä tzle processing for Toll activation) or binding to PRRs on the surface of immune tissues/cells (e.g., PGRP-LC binding for Imd activation in Drosophila). As shown in Table 3, transcripts of most Toll and Imd pathway proteins, such as Relish, Pelle, Cactus and Toll receptor, were influenced by parasitization. These include Toll proteins (GAJS01021639, log 2 Ratio PF/CF: -1.7) and Pelle (GAJS01022106, log 2 Ratio PF/CF: -1.5, log 2 Ratio PH/CH: -1.2) were down-regulated by parasitization (Table 3), which differs from P. xylostella after parasitization by D. semiclausum [12]. Overproduction of effector proteins, particularly anti-microbial peptides (AMPs), that immobilize pathogens, block their proliferation, or directly kill them is a hallmark of insect immunity [38,52]. Consistent with this notion, we have detected some AMPs and lysozyme. Most of them were up-regulated by parasitization (Table 3). However, it's worth pointing out that some immune response, like AMP genes, may be induced only by a puncture. Hence, some of the presented results may not be related to parasitism.
We also recorded changes in other proteins that influence immune responses in other moths such as tyrosine hydroxylase and dopa decarboxylase ( Table 3). The general finding is that expression of tyrosine hydroxylase and dopa decarboxylase is significantly induced following infection [12,53], while we found that tyrosine hydroxylase was up-regulated and dopa decarboxylase was down-regulated (Table 3).

Transcripts related to development and metabolism
Our data indicates that parasitism leads to up-regulation of genes associated with JH binding or degradation. JHBP (GAJS01017916, log 2 Ratio PF/CF: 1.6), JHE (GAJS01005072, log 2 Ratio PF/CF: 1.6) and JHEH (GAJS01070607, log 2 Ratio PF/CF: 1.6) ( Table 4). JHEH transcript levels were downregulated more than 2-fold in P. xylostella after parasitization by D. semiclausum [12]. Generally, JH is maintained at high levels during parasitoid larval development [5,[54][55][56]. Our findings run otherwise, with increases in JHE and JHEH transcript levels. This may be another example of the wide variation in molecular details of insect host-parasitoid relationships.
Parasitization of P. xylostella by D. semiclausum leads to downregulation of genes associated with ecdysteroid activities [12]. Our data support this view as the transcript level of ecdysteroid regulated protein (GAJS01010696, log2 Ratio PF/CF: -6.8) was down-regulated in parasitized larvae (Table 4). In general, parasitization leads to reductions in ecdysteroid titres and decreases in ecdysteroid regulated proteins [43,54,55,57]. We found that expression of methionine-rich storage protein (MRSP), a diapause-associated protein [58], was up-regulated (GAJS01008390, log 2 Ratio PF/CF: 6.4; log 2 Ratio PH/CH: 7.3) in parasitized larvae (Table 4). It indicates the crosstalk between regulatory pathways in insect diapause and in parasitoidregulated host development. We also found that the expression of transcripts encoding a number of G-protein coupled receptors (GPCRs) was affected by parasitization. For example, allatostatin  receptor (GAJS01001502, log 2 Ratio PF/CF: 4.7; log 2 Ratio PH/ CH: 2.1), neuropeptide A20 (GAJS01016991, log 2 Ratio PF/CF: 2.2), neuropeptide B3 (GAJS01001291, log 2 Ratio PF/CF: -3.0) and Methuselah (Mth) (GAJS01005512, log 2 Ratio PF/CF: -3.4; log 2 Ratio PH/CH: -1.8) were up-or down-regulated by parasitization (Table 4). We found that transcripts for the allatostatin receptor were highly expressed and it was reported that activation of allatostatin A-expressing neuron promotes food aversion and/or exerts an inhibitory influence on the motivation to feed in adult Drosophila [59]. Reduced food consumption was also found in parasitized larvae in our system [17]. There might be relationships among some GPCRs. Mth is a class B secretin-like GPCR and down-regulation of Mth increases the life span of D. melanogaster [60][61][62][63]. We infer that decreased Mth transcription levels may help elongate the lifespan of parasitized larvae, as seen elsewhere [1].

Validation of candidate genes
To validate our RNA-seq data, we performed qRT-PCR on 14 selected immune-and development-related genes, including ten anti-microbial peptides, MRSP, neuropeptide A1, serine protease inhibitor 7 and PAP 3 ( Figure 5, 6). The sequences of the primers used are given in Table 5. Our results are consistent with the RNA-seq profiles showing similar trends in up-or down-regulation (or steady-state) of host genes with little difference (Figure 5, 6). For example, based on RNA-seq analysis, attacin1 (q 6-fold), defensin (q 2-fold), MRSP (q 6-fold) and PAP 3 (q 2-fold) were upregulated in fatbody (Table 3, 4) and showed parallel changes in our qRT-PCR analysis ( Figure 5, 6). These data are based on pools of RNA from various time points. To investigate the expression of host genes at different periods after parasitization, we isolated RNA from 4 th instar larvae at selected time points after parasitization, and analyzed two associated genes, MRSP and attacin 1 (Figure 7). The expression of MRSP increased with time through 48 h following parasitization in fatbody; whereas transcripts of MRSP increased to a maximum at 24 h and then returned to the basal level by 48 h in hemocytes. Fatbody expression of attacin 1 peaked at 6 h p.p. and at 12 h in hemocytes.
We detected a range of CchBV genes in parasitized larvae. Based on our analysis, 19 unique sequences were identified from six PDV gene families including ankyrin, CrV1 protein, cystatin, EP, lectin and PTPs (Table 6). Besides these, eight other hypothetical proteins with unknown function were also found in parasitized larvae, and showed more than 47% similarity with some parts of BV reference genomes. The presence of multiple sequences for each PDV gene family is expected. We identified 2 segments with ankyrin domain, which are commonly shared by BV and IV, and 2 CrV1 transcripts, which showed high similarity (. 90%) with Cotesia sesamiae CrV1 protein [68] (Table 6). These two transcripts may belong to one gene. We detected two EP-like proteins. The EP genes encode secreted, glycosylated proteins that are expressed within 30 min in host and accumulate to comprise over 10% of total hemolymph proteins by 24 h p.p. [69], which could induce significant reduction in total hemocyte numbers and suppress host immune response presumably by its hemolytic activity during parasitization [70]. Three PTPs transcripts, which have 33 members as the largest family in CvBV [15], were also found in the parasitized larvae. Previous studies suggested that some but not all PTPs function as a phagocytic inhibitor or apoptosis inducer, which plays an important role in suppressing insect immune cell [71,72]. We found one lectin and one cystatin transcript in CchBV. All these proteins were also found in other BVs, like CvBV, CcBV and CrBV.
Transcription levels differed among different members and tissues with each gene family. Ankyrin 6, 2 and PTP1, 2, 3 had the lowest transcription levels relative to other CchBV genes in PCs-FB and PCs-HC (Figure 8). Ankyrin and PTP1 were mostly detected in PCs-FB. However, CrV1, cystatin 2 protein and lectin were mostly transcribed in PCs-HC and showed a higher transcription levels. EP1-like and EP2-like gene showed a parallel expression in both of infected tissues. Among the Unknown proteins, the unknown protein 2 and 4 had the highest transcription levels and expressed much more in PCs-HC ( Figure 8). The reasons why these genes expressed differential in different infected tissues need more investigation.

Conclusions
In summary, we investigated the global gene transcription profiles of fatbody and hemocytes of C. suppressalis in response to parasitization by C. chilonis. The results showed that the abundances of relatively low proportion of C. suppressalis transcripts were differentially expressed after parasitization by C. chilonis. And most of these affected genes have predicted roles in immunity, development, or detoxification. At the tissue level, our results indicated that the expressions of fat body genes were changed much more than those of hemocytes. Since we pooled the samples of four time point post-infection together, it is possible that we would loss the opportunity to assess patterns in transcriptome activity as function of sample time. Our results provide evidence for expression of 18 CchBV transcripts expressed in the host. The expression levels of these PDV genes were different at two tissues. It is also possible that CchBV gene products affect the growth and immune state of host through interactions at the protein level, such as viral proteins interacting with specific host proteins and epigenetic regulation. In addition, these viruses may also produce small non-coding RNAs that modulate host gene transcription or microRNA of host differentially expressed in response to parasitization [73]. The transcriptome data obtained in this study provides a basis for future research in this under-explored hostparasitoid interaction. Future functional studies on the identified immune-, development-and detoxification-related genes could lay the foundation for identifying hot-spots for host-parasitoid interaction, which could contribute to develop new strategies to optimize use of parasitoids for C. suppressalis control. Figure S1 Histogram presentation of Gene Ontology classification. The results are summarized in three main categories: biological process, cellular component and molecular function. The right y-axis indicates the number of genes in a category. The left y-axis indicates the percentage of a specific category of genes in that main category. The main and specific categories are indicated on the x-axis. (TIF) Figure S2 Histogram presentation of clusters of orthologous groups (COG) classification. All putative proteins were aligned to the COG database and can be classified functionally into at least 25 molecular families. (TIF)