Digital gene expression analysis of Helicoverpa armigera in the early stage of infection with Helicoverpa armigera nucleopolyhedrovirus

Accepted Manuscript Digital Gene Expression analysis of Helicoverpa armigera in the early stage of infection with Helicoverpa armigera nucleopolyhedrovirus Zhen Li, Zhenqiang Lu, Xiu Wang, Songdou Zhang, Qingwen Zhang, Xiaoxia Liu PII: DOI: Reference:

S0022-2011(15)30007-0 http://dx.doi.org/10.1016/j.jip.2015.08.008 YJIPA 6724

To appear in:

Journal of Invertebrate Pathology

Received Date: Revised Date: Accepted Date:

8 May 2015 12 August 2015 18 August 2015

Please cite this article as: Li, Z., Lu, Z., Wang, X., Zhang, S., Zhang, Q., Liu, X., Digital Gene Expression analysis of Helicoverpa armigera in the early stage of infection with Helicoverpa armigera nucleopolyhedrovirus, Journal of Invertebrate Pathology (2015), doi: http://dx.doi.org/10.1016/j.jip.2015.08.008

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

For Consideration: Journal of Invertebrate Pathology -------------------------------------------------------------------------------------------------------

Digital Gene Expression analysis of Helicoverpa armigera in the early stage of infection with Helicoverpa armigera nucleopolyhedrovirus

Zhen Li a,#, Zhenqiang Lu a,#, Xiu Wang a, Songdou Zhang a, Qingwen Zhang a, Xiaoxia Liu a,* a

Department of Entomology, China Agricultural University, Beijing 100193, China

------------------------------------------------------------------------------------------------------#

: These authors contributed equally to this study

* Corresponding Author: Dr. Xiaoxia Liu Department of Entomology China Agricultural University Yuanmingyuan West Road Beijing 100193, China Phone: 86-10-62733946 Fax: 86-10-62733946 Email: [email protected]

1

ABSTRACT Helicoverpa armigera single nucleocapsid nucleopolyhedrovirus (HearNPV) is an obligatory and lethal parasite of the cotton bollworm and has been extensively used in China for the control of this notorious pest. Digital gene expression (DGE) analysis was adopted for an overall comparison of transcriptome profiling between HearNPV-infected and control healthy H. armigera larvae during an early stage post-inoculation. A total of 908 differentially expressed genes (DEGs) were identified, of which 136 were up-regulated and 597 were down-regulated. GO category and KEGG pathway analysis demonstrated that the identified DEGs involved in ribosome biogenesis, aminoacyl-tRNA biosynthesis, protein processing in endoplasmic reticulum, biosynthesis of valine, leucine, isoleucine and the spliceosome were significantly down-regulated, whereas genes involved in pancreatic secretion, protein digestion and absorption and salivary secretion showed obviously up-regulated transcription. The DEGs were verified by quantitative real-time PCR, and genes that participated in defensive response, nutritional digestion and developmental regulation exhibited specific expression patterns in a continuous time-course assessment. These results provide basic data for future research on the molecular mechanism of HearNPV infection and the interactions between lepidopteran hosts and their specific NPV parasites.

Keywords: cotton bollworm, early response, differentially expressed genes, quantitative real-time PCR, gene network

2

1. Introduction

A number of nucleopolyhedroviruses (NPV) belonging to family Baculoviridae are lethal pathogens of phytophagous pests of Lepidoptera, Hymenoptera and Diptera in fields (Chen et al., 2001; Moscardi, 1999). They typically cause the death of the host larvae 5 to 7 days post-infection. Additionally, NPV are normally obligatory for specific pests and can be rapidly transmitted through aerial flow and propagation carriage; these properties have been commercialized and widely used as environment-friendly substitutes for chemical pesticides (Grzywacz et al., 2000). Helicoverpa armigera single nucleocapsid nucleopolyhedrovirus (HearNPV) belongs to the genus Alphabaculovirus and is an obligatorily parasite of the cotton bollworm. HearNPV has been extensively used in China for the control of this notorious pest in cotton and vegetable fields (Chen et al., 2001; Georgievska et al., 2010). HearNPV contains a single nucleocapsid packaged into a virion, which is then occluded in a polyhedron. After being fed on by lepidopteran larvae, the NPV polyhedral inclusion bodies are digested by the alkaline midgut fluid. The virions are subsequently released into the epithelial cells of the midgut and then spread throughout the insect body (Chen et al., 2001; Etebari et al., 2007). Finally, the diseased larvae are totally liquefied and left covered with only a fragile cuticle, and the corpse serves as a virus donor for further transmission in the field. Currently, studies on the pathogenesis of NPV disease are mainly focused on the investigation of auxiliary genes encoded by NPV and their facilitating functions in infection, host body liquefaction and host behavior manipulation (Hawtin et al., 1997; Hoover et al., 2011; Kamita et al., 2005; Lanier et al., 1996). In contrast, studies on the early response are relatively limited, especially studies investigating the differentially expressed genes in insects upon NPV infection and their potential involvement in the interactive mechanisms. The development of high-throughput deep sequencing technologies has had a dramatic impact on the analysis of differential gene expression among different samples (Morozova and Marra, 2008) and has accelerated the broad understanding of host responses to pathogen infections at a molecular level (Gao et al., 2014; Hou et al., 2014; Yue et al., 2015). Comparative expressed-sequence tag analysis and genome-wide analysis of gene expression in silkworm cells 3

infected with Bombyx mori nucleopolyhedrovirus (BmNPV) were conducted to estimate the host response to NPV infection (Okano et al., 2001; Sagisaka et al., 2010). A similar study was also performed

in

Spodoptera

frugiperda

cells

infected

by

Autographa

californica

nucleopolyhedrovirus (Nobiron et al., 2003). Recently, differential mRNA expression was analyzed in Helicoverpa zea cells infected with a baculovirus of H. armigera (Nguyen et al., 2013). These studies all assessed alterations in the expression of genes in insect cell lines infected with NPV. However, insects are complex organism, and their responses to pathogen infection usually result from integrated reactions from different tissues and systems in the whole insect body and not only from specific cells. Digital gene expression (DGE) is one of the recently advanced sequencing technologies. DGE can screen and identify millions of differentially expressed genes (DEGs) in samples without prior annotations, and thus has been widely adopted for transcriptome profiling studies (Bennett et al., 2005). In the current study, DGE with the Illumina Genome Analyzer platform was adopted for a systematical assessment between HearNPV-infected and control healthy H. armigera larvae at 24 hpi (hours post-inoculation) to obtain an overall transcriptome profile of whole H. armigera larvae responding to early infection with HearNPV. The up- and down-regulated DEGs were seriously determined, their Gene Ontology (GO) and KEGG pathways enrichments were analyzed; qRT-PCR was conducted to detect and verify the gene expression levels of the identified DEGs, and potential networks of proteins encoded by the identified DEGs were predicted. A large number of transcripts were found to be down-regulated upon NPV-infection, and genes involved in the defensive response, nutritional digestion, and growth and development regulation showed obvious expression patterns in the time-course analysis. The present study may supply information to help elucidate the response of H. armigera larvae to infection with HearNPV and the interactive mechanism between NPV and their hosts at a transcriptional level.

2. Materials and methods

2.1. Insect and virus Larvae of the cotton bollworm Helicoverpa armigera Hübner (Lepidoptera: Noctuidae) were 4

originally collected from cotton fields in Hebei Province, China, and were reared under laboratory conditions (27±1 °C, 80±10% RH, 14 h L: 10 h D) with an artificial diet (Fan, 2003) for at least thirty generations. For larvae preparation for the following tests, healthy late third instars were selected and individually collected in finger-shaped glass tubes (5.5 cm in length × 2.2 cm in diameter) without food and starved for 8 h until molting. Raw powder of HearNPV occlusion bodies (OBs; 5 × 1011 OBs g-1) was purchased from the Henan Jiyuan Baiyun Industry Co., Ltd (Henan, China) and stored at 4 °C for later use.

2.2 Toxicity test A total of 0.02 g of HearNPV raw powder was thoroughly suspended in 1 ml of sterilized ddH2O by vortexing for 5 min. The resulting product was used as a virus stock with a HearNPV concentration of 1 × 1010 OBs ml-1. For the toxicity test, 1 × 10 to 1 × 109 OBs ml-1 of the HearNPV suspensions was prepared by 10-fold serial dilutions of the stock. For oral inoculation of the virus, 10 µl of the HearNPV suspension was pipetted onto one piece of artificial diet jelly (8 mm in length × 8 mm in width × 1 mm in thickness), resulting in 1 × 10-1 to 1 × 107 OBs of HearNPV per diet piece. Then, one piece of the diet with virus was placed into a finger-shaped glass tube together with one newly molted fourth instar after starvation pretreatment. A diet piece with the addition of 10 µl of sterilized ddH2O was supplied as a control (0 OBs treatment). Subsequently, only the larvae that ate up the virus-contaminated or control diet were selected, replenished with a sufficient fresh diet without virus, reared individually under normal conditions and monitored at 8 h intervals for mortality. The survey lasted until the tested larvae died or pupated. The lethal dose value (LD50) was determined by probit analysis between the probability of death and the dose of HearNPV using SPSS software (IBM SPSS statistics 20).

2.3. Confirmation of HearNPV infection To track and verify the infection of H. armigera with HearNPV, the larvae were inoculated with 1 × 105 OBs of HearNPV. Subsequently, the head (H), midgut (MG), ventral nerve cord (VNG), fat body (FB) and Malpigian tube (MT) of the HearNPV-inoculated and control larvae were dissected in cold PBS (pH 7.8) at 24, 48, 72 and 96 hpi, and the transcriptional level of 5

polyhedrin (polh, NP_075070), a very-late expressed gene of NPV that encodes the polyhedrin for the occlusion body (Iatrou et al., 1985) was quantified using qRT-PCR. The primers are listed in Table S1 and the specific qRT-PCR procedure is described below.

2.4. Total RNA extract HearNPV-inoculated or control larvae were randomly collected and snap-frozen in liquid nitrogen. Five pooled frozen larvae were ground into powder in liquid nitrogen with mortars and pestles; then, the powdered aliquots were transferred into 1.5 ml RNase-Free Eppendorf tubes. The total RNA from the samples was isolated using the Qiagen RNeasy Mini Kit (Qiagen, Valencia, CA, USA). The contaminated genomic DNA was digested with RNase-Free DNase (Qiagen). The RNA samples were qualified and quantified using the Nanodrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) and Agilent bio-analyzer (Agilent, Santa Clara, CA, USA).

2.5. DGE library construction and sequencing DGE libraries (HearNPV-inoculated and control) were constructed using the Illumina Gene Expression Sample Prep Kit as described in our previous study (Yang et al., 2014). Briefly, mRNA was enriched from 6 µg of total RNA extracted from the control and NPV-inoculated larvae that ate the virus diet piece using oligo(dT) magnetic beads at 24 hpi. The first strand cDNA was synthesized by random hexamer-primers from the template of mRNA short fragments (approximately 200 bp) that were interrupted from the enriched mRNA using fragmentation buffer. After the second strand of cDNA was synthesized by DNA polymerase I, purified, end repaired and adenine was added, sequencing adaptors were ligated to the fragments. Finally, the libraries were constructed and used for sequencing after fragment purification by agarose gel electrophoresis and enrichment by PCR amplification. DGE library products were sequenced via the Illumina HiSeqTM 2000.

2.6. Determination of DEGs Raw reads were obtained after sequencing, and clean reads were generated by filtering and

6

removing dirty raw reads including adaptors, reads with unknown bases and low quality reads. The resulting clean reads were aligned to reference sequences (B. mori genome and H. armigera ESTs) using SOAP aligner/-soap2 (Li et al., 2009) with mismatches of no more than 2 bases allowed in the alignment. Then, the sequencing data were measured through analysis of sequence saturation and read distributions on the reference genome/ESTs. Statistical analysis was performed to summarize the number of clean reads that aligned to the reference genome/ESTs, and the gene expression of transcripts was calculated by the numbers of reads mapped to the H. armigera ESTs using the RPKM (Reads Per kb Million Reads) method (Mortazavi et al., 2008). Then, the differentially expressed genes between HearNPV-inoculated and control larvae were determined according to the algorithm developed by Audic and Claverie (1997) as previously described (Yang et al., 2014). An FDR (False Discovery Rate) ≤ 0.001 and an absolute value of differential expression level ≥ 2-fold were used as criteria to identify DEGs.

2.7. Gene ontology and pathway functional enrichment analysis of DEGs Functional annotation of DEGs was conducted by BLASTx analysis against the NCBI non-redundant protein database with a cut-off E-value of 10-5 (Yang et al., 2014). Then, the BLASTx

results

were

mapped

to

Gene

Ontology

(GO)

terms

in

the

database

(http://www.geneontology.org/), and the significantly enriched GO terms in the DEGs were calculated using the hypergeometric test with a corrected P value ≤ 0.05 as the threshold. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed against the public pathway-related database (http://www.genome.jp/kegg ) (Kanehisa et al., 2008). Pathways containing DEGs with P values ≤ 0.05 were determined to be significantly enriched metabolic or signal transduction pathways.

2.8. Quantitative real-time PCR analysis To validate the DEGs from Illumina Sequencing, twenty DEGs with∣log2 Ratio∣≥2 were randomly selected and verified with qRT-PCR. The sequences of the selected transcripts were retrieved from the EST library of H. armigera from the NCBI database. Eleven genes from the

7

above twenty DEGs of interest were quantified with a time course analysis by qRT-PCR and compared between the HearNPV-inoculated and control larvae at 0, 6, 12, 24, 48, 72, and 96 hpi. RNA extracted from a mix of at least 5 larvae was regarded as one biological replicate, and three independent biological replicaties were adopted for each treatment. Three technical repetitions were also performed for each treatment at each tested time point. Gene-specific primers for all of the tested genes were designed using the online software Primer 3 (http://www.simgene.com/Primer3) and Primer BLAST in the NCBI tools; the primers are listed in Table S1. An amplification system (20 µl) containing 10 µl of SYBR Premix Ex TaqTM (Takara, Otsu, Shiga, Japan), 4 µl of each specific primer (1 µM) and 1 µl of cDNA (10-fold dilution) was prepared and tested on the CFX ConnectTM Real-Time PCR System (Bio-Rad, Hercules, CA, USA) using the following program: 95 °C for 30 s for pre-denaturation, 40 cycles of 95 °C for 5 s and 60°C for 30 s for PCR amplification and signal collection, and 95 °C for 15 s, 60°C for 60 s and 95 °C for 15 s for the dissociation stage. H. armigera genes including elongation factor 1α (EF1α) and ribosomal protein L13 (RPL13) were selected as reference genes after stability evaluation. The geometric mean calculated from the CT values of these two genes was used as a normalization factor for gene expression analysis using the ΔΔCT method. The significance of the expression difference between the HearNPV-inoculated and control larvae was statistically analyzed using the independent-sample T test with SPSS Statistics 17.0 (SPSS Inc., Chicago, IL, USA).

2.9 Gene network prediction The interactions among the identified DEGs (∣log2 Ratio∣≥2) were analyzed using the online database search tool GeneMANIA (http://www.genemania.org). The interactions among the up-regulated genes and the down-regulated genes were respectively predicted. The predicted gene networks were generated with Cytoscape 2.8.3.

3. Results

8

3.1. Selection of the HearNPV dose The mortality of H. armigera followed an S-shaped curve. The regression analysis illustrated that the LD50 of HearNPV for H. armigera larvae was 12,046.93 OBs, with 95% confidence limits from 6046.93 OBs to 20,950.35 OBs. To investigate the early response of the host H. armigera at the transcriptional level, 105 OBs of HearNPV (which could induce approximately 51.3% mortality of fourth instars and 82.9% mortality in total) was selected for oral inoculation and DGE analysis.

3.2. Confirmation of HearNPV infection The qRT-PCR demonstrated that the expression of NPV-encoded polh increased with the time post-inoculation and dramatically increased from 48 hpi in all tested tissues. HearNPV showed the highest expression level in the fat body at 96 hpi (Fig. 1).

3.3. Profiles of the DGE libraries To elucidate the early gene transcription response in H. armigera, two DGE libraries were constructed from the uninfected and HearNPV-inoculated H. armigera larvae at 24 hpi using Illumina sequencing technology. Table.1 demonstrated the statistics of the DGE reads. A total of 4,842,523 and 4,765,688 raw reads were generated in the control and NPV-inoculated DGE libraries, respectively. After removing reads with low quality, containing N’s and only adaptors, 4,817,921 and 4,736,444 clean reads accounting for more than 99.3% of the raw reads were filtered from the two libraries. In the control and HearSNPV-inoculated libraries, 26.18% and 43.28% of the clean reads corresponding to 22.31% and 37.24% of the unique reads could be mapped to genes in the H. armigera ESTs, while only 2.99% and 5.36% of the clean reads corresponding to 2.53% and 4.47% of the unique reads were mapped to genes in the B. mori genome. Therefore, the H. armigera ESTs were used as reference data for the following DEG analysis. To estimate the sequencing depth, the sequencing saturation was analyzed. The results showed that the percentage of identified genes increased along with the number of clean reads, and the

9

number of identified genes approached saturation when the number of total reads reached approximately 2M (Fig. S2). Our sequencing depth was approximately 5M in both DGE libraries (Fig. S2), which was sufficient transcript coverage.

3.4. Analysis of differentially expressed genes RPKM values were calculated to estimate the abundance of unique reads in the two libraries. In the control library, 75.8% of the reads had RPKM values more than 50, 23.4% were between 10 and 50, and only 3.8% were below 10. However, in the HearNPV-inoculated library only 43.7% of reads had RPKM values more than 50, 43.6% were between 10 and 50, and 12.5% were less than 10 (Table. S2). Comparatively, the distributions of gene abundance were different between the two libraries. The number of genes expressed at a high level (above 50) in the control library was increased compared to the HearNPV-inoculated library, but more genes with lower expression levels were detected in the HearNPV-inoculated library. Identification of DEGs prior to the widespread propagation of HearNPV in the host may provide important clues to the interactions between the host and virus, including the host response against NPV infection and the mechanism of NPV infection. DEGs were determined between the control and NPV-inoculated libraries using an FDR ≤ 0.001 and∣log2 Ratio∣≥1 as a threshold. The analysis showed that a total of 908 genes showed altered gene expression levels, with 190 up-regulated and 718 down-regulated (Fig. S3).

3.5. GO enrichment analysis To clarify the functions of the DEGs, all 908 identified DEGs were mapped to terms in the GO database. The analysis results showed that 136 up-regulated genes and 597 down-regulated genes were annotated and categorized into three ontologies: molecular function, cellular component, and biological process. Among the 136 functions of the known up-regulated genes, 47 genes were classified into four subgroups of biological process [i.e., metabolic process (34), cellular process (2), biological regulation (5) and defense response (6)], 68 genes were categorized into four subgroups of molecular function [catalytic activity (39), binding (12), enzyme regulator activity

10

(13) and transporter activity (4)] and 21 genes were sorted to cellular component (Table. 2). Among the 597 functions of the known down-regulated genes, 191 genes were classified into five subgroups [metabolic process (56), cellular process (60), biological regulation (46), translation (14), and defense response (15)], 237 genes were categorized into four subgroups [catalytic activity (40), binding (157), enzyme regulator activity (37) and transporter activity (3)] and 169 genes were sorted to cellular components (Table. 2). Comparatively, the translation category was only detected in the down-regulated DEGs. The terms “catalytic activity” and “metabolic process” in the up-regulated DEGs and “binding” among the down-regulated DEGs were dominant. The sequences of 28 up-regulated and 62 down-regulated DEGs with∣log2 Ratio∣≥2 and FDR ≤ 0.001were retrieved in the NCBI nucleotide database according to their accession numbers and were listed in Table. 3 and Table. 4. Genes involved in “catalytic activity” and “defense response” represented the majority of the up-regulated DEGs, while genes termed in “biological regulation” represented the majority of the down-regulated DEGs. The homologs of DEGs were obtained after performing a standard nucleotide blast using the basic local alignment search tool of NCBI, and their corresponding accession numbers, highest alignment scores and best expected values (E values) were also described.

3.6. KEGG pathway analysis The KEGG pathway database, which records networks of molecular interactions in the cell and their variants that are specific to particular organisms, was used to increase our understanding of the biological functions of the identified DEGs in the H. armigera larvae before widespread propagation of the inoculated HearNPV. A total of 217 genes consisting of 55 up-regulated genes and 162 down-regulated ones were mapped onto 15 significant (P < 0.05) KEGG pathways (Table. 5). The down-regulated DEGs were significantly mapped to all 15 pathways, while the up-regulated DEGs were only significantly enriched in 7 KEGG pathways. The down-regulated DEGs showed more involvement in the 5 pathways of protein biosynthesis and translation (ribosome biogenesis, aminoacyl-tRNA biosynthesis, protein processing in endoplasmic reticulum, biosynthesis of valine, leucine, isoleucine and spliceosome), 2 pathways of response to infection (p53 signaling pathway and antigen processing and presentation), 1 pathway of pyrimidine 11

metabolism, and 1 pathway involving neuroactive ligand-receptor interactions. The up-regulated DEGs showed more participation in the 3 pathways of nutrient digestion (pancreatic secretion, protein digestion and absorption and salivary secretion).

3.7. Validation of differentially expressed genes The difference in transcriptional levels between the control and NPV-inoculated samples at 24 hpi was determined by qRT-PCR and compared with the corresponding expression profiles obtained from DEG analysis. The validation of all twenty tested genes using qRT-PCR showed concordant directions of differential gene expression determined by DGE analysis, with only some fold-change discrepancies (Fig. 2). To obtain a better understanding of the DGE results, eleven genes of interest involved in the defense response, nutritional digestion and developmental regulation were selected for time-course analysis at 0, 6, 12，24, 48, 72, and 96 hpi using qRT-PCR detection. Two genes [stress-induced phosphoprotein 1(stip1) and heat shock protein 70 (hsp70)] that were involved in the defense response showed induced expression levels at most of the time points (Fig. 3A). However, the transcriptional level of serine protease 52, which is also an important protein in immunity, was significantly decreased throughout the investigation (Fig. 3A). The genes participating in nutrient metabolism, including glucose oxidase-like enzyme (gox), fructosidase, carboxyl/choline esterase (cce) (Fig. 3B), trypsin carboxypeptidase precursor and neutral lipase (Fig. 3C) exhibited induced expressional patterns. The ecdysteroid-regulated protein and juvenile hormone diol kinase genes, which are related to the regulation of growth and development, showed induced transcription patterns after HearNPV infection (Fig. 3D).

3.8 DEG networks The retrieval of gene interactions with GeneMANIA found potential networks among 3 up-regulated and 26 down-regulated DEGs (∣log2 Ratio∣≥2) (Fig. 4). The predicted genetic interactions of the up-regulated DEGs happened among hydroxyacid oxidase 1-like (HAO1), peptidoglycan recognition protein B (PGRPb) and sepiapterin reductase (SPR) (Fig. 4A). The 12

gene networks of the down-regulated DEGs were constructed among genes primarily participating in: 1) the process of DNA replication and transcription, RNA translation and protein conformation (the

majority

of

DEGs

in

the

network);

2)

detoxification

and

immune

[i.e.,

choline/ethanolaminephosphotransferase 1-like (CEPT1), protein Mpv17-like (MPV17), and stress-induced phosphoprotein 1 (STIP1)]; 3) cytoskeletal construction [talin-2-like (TLN2)] (Fig. 4B).

4. Discussion

HearNPV is one of the most widely used commercialized NPV products for the control of the devastating pest H. armigera in cotton and vegetable fields (Georgievska et al., 2010; Grzywacz et al., 2000). The pathogenicity of HearNPV declines when H. armigera larvae grow larger or the inoculation dose decreases. In our study, the newly molted 4th instars showed moderate susceptibility to HearNPV infection and were determined to be the best developmental stage for the study of the interaction between HearNPV and H. armigera larvae. To ensure NPV infection and select an observable time for the detection of the host response, 105 OBs of HearNPV (which could kill most of the host larvae within 5 to 7 days post-inoculation) was selected as the inoculum. Typically, host larvae show symptoms from the infection and NPV multiplication from the second to fourth day post-infection, resulting in a reduction in food intake and molting time (Riegel and Slavicek, 1997). In the present study, a dramatic increase in polh transcripts encoded by HearNPV was found in different tissues of the infected larvae from 48 hpi and reached the highest expression level in the fat body at 72 hpi. Previous analysis of mRNA expression of H. zea insect cells infected with HearNPV demonstrated that 24 hpi was a boundary point when the expression of genes encoded by NPV shifted from a small number to the majority (Nguyen et al., 2013). Moreover, gene transcription and protein synthesis in the host larvae would be shut off by the increased multiplication of NPV in the host (Nobiron et al., 2003; Ooi and Miller, 1988; Popham, et al., 2010). Therefore, 24 hpi was selected as the sampling time for the following DGE analysis to understand the early response of H. armigera larvae to HearNPV infection before the 13

wide-spread propagation of the virus. Because only a small number of clean reads could be mapped to the B. mori genome, the H. armigera EST database was chosen as the reference for the DGE analysis. The large number of DEGs and noticeably higher numbers of down-regulated DEGs (718) compared to the up-regulated DEGs (190) determined in current study were in accordance with a previous study on DEG analysis in H. zea insect cells infected with HearNPV in which more number of DEGs were found at 24 hpi, with 305 up-regulated vs 19673 down-regulated genes (Nguyen et al., 2013). Autographa californica multiple nucleopolyhedrovirus (AcMNPV) infection in Trichoplusia ni also caused a dramatic down-regulation of unigenes (3039) compared with the 32 up-regulated unigenes detected at 24 hpi (Chen et al., 2014). However, only 52 DEGs were determined in the BmNPV-infected B. mori cell lines at 12 hpi, and the number of up-regulated genes (35) was higher than the number of down-regulated genes (17) (Sagisaka et al., 2010). At the translation level, 18 significantly up-regulated vs 12 down-regulated proteins were detected at 24 hpi in Hz-AM1 cells (from Heliothis virescens) infected with AcMNPV (Popham et al., 2010). In larvae of the light brown apple moth Epiphyas postvittana, notably more up-regulated genes (84) were also determined than down-regulated genes (18) in response to Epiphyas postvittana nucleopolyhedrovirus (EppoNPV) infection on the fifth day post-inoculation (dpi) (Gatehouse et al., 2009). These results suggested that the interactions between variant insects and their corresponding NPV might be divergent; however, most host genes showed a down-regulated expression pattern over the time course of infection, with only a small number of up-regulated genes detected. The down-regulated expression of the majority of H. armigera genes 24 h post-HearNPV infection might lead to the collapse of the host defensive response, probably due to NPV manipulation, and would benefit successful multiplication and colonization of HearNPV in H. armigera. GO and KEGG analysis were performed for the investigation of DEG enrichment. In the reported systems of NPV-insect interactions, genes with roles in the stress response, apoptosis, immunity, transcription and translation regulation, protein trafficking, signal transduction, energy generation, endocrine modulation and cytoskeleton rearrangement were primarily enriched and clustered (Bao et al., 2009; Chen et al., 2014; Gatehouse et al., 2009; Nguyen et al., 2013; Popham 14

et al., 2010; Sagisaka et al., 2010; Salem et al., 2011). In this study, genes termed in “catalytic activity” and “metabolic process” were primarily induced by infection with HearNPV, which might facilitate NPV infection and multiplication due to increased nutrition and energy supply through digestion from host storage. Genes involved in sugar transportation and carbon metabolism were also found to be up-regulated in the interactions of BmNPV-B. mori and AcMNPV-S. frugiperda (Choi et al., 2012; Nguyen et al., 2013; Sagisaka et al., 2010); genes that encoded proteins associated with mitochondria and energy/metabolism, including NADH dehydrogenase and ATP synthase, were also induced in the interaction of AcMNPV-T. ni (Chen et al., 2014). No NPV-encoded genes have been linked to energy metabolism (Chen et al., 2001; Nguyen et al., 2013). Similarly, the baculovirus manipulation of the host’s metabolic machinery has been proposed to be critical for their successful propagation in hosts (Monteiro et al., 2012). “Defensive response” is one of the most considerable reactions in the host upon pathogen infection, but most of the DEGs termed in this category were found to be down-regulated in this study; the exceptions were of genes encoding the D-type cecropins-like, peptidoglycan recognition protein B, attacin and gloverin (Table. 3). The immune response gene cecropin was also found to be up-regulated in HearNPV-H. zea, AcMNPV-S. frugiperda and other DNA virus-insect systems (Choi et al., 2012; Nguyen et al., 2013). The expression level of gloverin (a lepidopteran-specific antibacterial peptide) in the midgut of B. mori could be strongly induced during the early stage of BmNPV infection (Bao et al., 2009). The 4.22-fold (log2 ratio) up-regulation of the gene encoding D-type cecropins-like and 2.65-fold (log2 ratio) up-regulation of gloverin detected in this DGE analysis suggested that these two genes probably served as important host defense reagents against baculovirus. KEGG analysis identified 15 significantly enriched pathways in the DEGs. Most of the down-regulated DEGs were associated with the ribosome biogenesis pathways. Ribosome biogenesis in eukaryotes takes place both in the cytoplasm and the nucleolus and is closely linked to other cellular activities, such as growth and division (Thomson et al., 2013). The Notch signaling pathway, which participates in the regulation of the process of cell division, was also consistently suppressed during early infection with HearNPV in the H. zea cell line (Nguyen et al., 2013). Hence, the interruption of cell division via the suppression of ribosome biogenesis and 15

signal transduction might be one possible mechanism underlying the delayed growth of HearNPV-infected H. armigera. Moreover, DEGs in other significantly enriched pathways involved in protein biosynthesis and translation (i.e., aminoacyl-tRNA biosynthesis, protein processing in endoplasmic reticulum, biosynthesis of valine, leucine, isoleucine and spliceosome) were also found to be down-regulated in this study. Similarly, upon infection with EppoNPV, some protein production-related genes (i.e., translation initiation factor and elongation factor) were down-regulated at the transcriptional level in E. postvittana (Gatehouse et al., 2009). Therefore, interference with normal biosynthesis and cell growth and division in the H. armigera larvae probably represent important steps in the early infection of HearNPV. Pathways involved in the defensive reactions in the NPV-insect system have been a primary focus in previous studies. In the present analysis, two defense-related pathways [p53 signaling pathway (also named the p53 tumor suppressor pathway) and antigen processing and presentation pathway] were significantly enriched at 24 hpi, and the majority of genes in these two pathways were notably down-regulated. At an earlier infection phase (18 hpi) of HearNPV in H. zea cells, the apoptosis enhancer p53 gene was found to be up-regulated and the apoptosis inhibitors (including the IAP gene) were down-regulated, which triggered apoptosis in host cells for defense against virus infection (Nguyen et al., 2012). Their increased transcriptional levels were reduced concomitant with the infection (Chen et al., 2014), resulting in the collapse of host defense reactions and successful NPV infection. Manipulation of host defense reactions may be an important strategy for HearNPV and is critical for further colonization. DEGs in pathways involved with nutrient digestion in H. armigera (i.e., pancreatic secretion, salivary secretion and protein digestion and absorption) were highly

up-regulated.

Sugar

metabolism-related

genes,

including

genes

encoding

UDP-glucosyltransferase and glucosyl/glucuronosyl transferase, were demonstrated to be up-regulated in EppoNPV-infected E. postvittana (Gatehouse et al., 2009). Genes involved in protein degradation and processing were up-regulated in the fat body of susceptible B. mori strain than in resistant strain (Bao et al., 2010). Genes clustered in protein metabolism were reported to be up-regulated in H. zea cells upon HearNPV infection at 18 hpi (Nguyen et al., 2012). Therefore, we hypothesized that the virus could abrogate the defensive reactions during the preliminary stage of HearNPV infection, inhibit normal biosynthesis and accelerate nutrient digestion in the 16

H.armigera host. These effects would accommodate further multiplication and colonization of HearNPV with suitable and nutritional circumstances. QRT-PCR has been considered to be an important complements to DGE analysis for validation. In this study, all of the genes validated with qRT-PCR showed expression patterns that were similar to the DGE results, suggesting that DGE was a credible technique for the large-scale screening of differentially expressed candidate genes among samples. Compared with qRT-PCR, DGE analysis for limited replication may over or under estimate the fold changes of expression differences of some tested genes. In the current study, the fold changes of genes encoding neutral lipase, trypsin, juvenile hormone diol kinase, serine protease 52, and heat shock protein 90 were more intensive in the qRT-PCR than in the DGE analysis. Many reports have shown this type of divergences between DGE and qRT-PCR analysis (Gao et al., 2014; Yang et al., 2014; Yue et al., 2015). Therefore, qRT-PCR validation is indispensible for sophisticated results concerning specific genes of interest. The influences of NPV infection on gene expression in the host larvae change concomitant with the process of virus infection. Indeed, the majority of gene expression could be shut off from 24 hpi (Nguyen et al., 2013). Significantly more down-regulated genes were also detected in HearNPV-infected larvae at 24 hpi based on the DGE analysis, whereas in the 96 hour time-course analysis of gene expression in the current study, several genes in H. armigera were found to have an induced expression pattern that was concomitant with HearNPV infection; this result was probably due to NPV manipulation to aid virus infection and multiplication. Hsp70, which exhibited continuously induced transcription in the diseased larvae, has been widely studied in varieties of NPV-insect interactive systems (Breitenbach and Popham, 2013; Iwanaga et al., 2014). In addition to its stress response activity and chaperone role for protein folding, HSP70 was proposed to be a novel structural protein of BmNPV and provide a clue to the interaction between the host defensive response and NPV multiplication (Iwanaga et al., 2014). Thus, the continuously induced expression level of hsp70 in the later phase (48-96 hpi) of HearNPV infection might be critical for the multiplication and infection of HearNPV in H. armigera. The induced expression of genes that participate in nutrient digestion and metabolism was also detected in the time-course analysis. Previous researchers found that baculoviruses could successfully use the metabolic 17

machinery of its host for fast multiplication (Maynard et al., 2010; Monteiro et al., 2012). Therefore, the manipulated acceleration of food digestion and transportation in diseased H. armigera might accommodate the propagation and infection of HearNPV. Delayed molting and prolonged larval duration of HearNPV-infected H. armegera larvae were detected in our study. The baculovirus-encoded gene ecdysteroid uridine 5’-diphosphate (UDP)-glucosyltransferase (egt) has been revealed to play critical role in the decreas of active ecdysteroid levels in the hemolymph of AcMNPV-infected S. frugiperda and LdMNPV-infected Lymantria dispar. This decrease resulted in molting prevention and up-climbing behavior initiation of the infected host larvae (Hoover et al., 2011; O’Reilly et al., 1992). In addition to the effect on ecdysone metabolism, a juvenile hormone binding protein was also detected with up-regulated transcripts at 24 hpi in AcMNPV-infected S. frugiperda (Salem et al., 2011). In HearNPV-infected H. armigera, insect developmental regulated genes including ecdysteroid-regulated protein and juvenile hormone (JH) diol kinase (JHDK) showed abnormal up-regulated expression. The disorder in the hormone level s in the diseased larvae and the resulting delayed molting and prolonged instar duration would benefit fast infection and multiplication of NPV, as the molting of insects is energy consuming and can exclude pathogens from the insect body (Hoover et al., 2011). Gene network provides information for better understanding of DEG functions. In the current study, network of the up-regulated DEGs were constructed among HAO1, PGRPb and SPR. HAO1 and PGRPb involve in detoxification and immune, and SPR participates in diminishing the levels of dopamine, norepinephrine and serotonin in organism. As biogenic amine, dopamine, norepinephrine and serotonin serve as important neurotransmitter in insects and participate in regulation

of

variant

physiological

states

and

behaviors

(Monastirioti,

1999).

The

BmNPV-infected silkworms showed enhanced locomotory behavior on the 4 dpi (Kamita et al., 2005), HearNPV-infected cotton bollworm (Georgievska et al., 2010) and LdMNPV-infected European gypsy moth (Hoover et al., 2011) all exhibited climbing behavior in the later stage after NPV infection. In previous studies, two NPV-encoded genes (ptp and egt) were found to be critical to the manipulated behavior of NPV-infected lepidopteran larvae (Hoover et al., 2011; Kamita et al., 2005). From the gene network predicted in present study, biogenic amine might be another factor in behavior regulation of NPV-infected insects. Concomitant with collapse of the 18

defense reactions in HearNPV-infected H. armigera, the homeostasis level of neurotransmitter in H. armigera might be destroyed and result in the behavior alteration of HearNPV-infected larvae in the later phase of HearNPV-infection. Consistent with the GO and KEGG analysis, the down-regulated DEGs showed primarily networks in the biosynthesis and conformation of proteins. The down-regulated protein process probably also resulted in the breakdown of the defense line in the HearNPV-infected H. armigera through networks with DEGs involved in detoxification and immune response. In summary, 908 DEGs were identified between control and HearNPV-inoculated cotton bollworm larvae at 24 hpi via DGE analysis, and most of the DEGs were down-regulated. GO terms, KEGG pathway analysis, qRT-PCR verification and gene network prediction demonstrated that the identified DEGs involved in protein biosynthesis and translation were significantly down-regulated, the DEGs involved in nutrient digestion showed obviously up-regulated transcription, and networks of genes happened among many DEGs. We concluded that the majority of genes and metabolic levels would be gradually shut off in H. armigera larvae upon HearNPV infection and NPV-manipulated alterations in nutrient digestion, defensive reactions and developmental processes of the host insects might be a critical strategy for virus multiplication and infection. The present study provides an overview of the gene transcriptional profile in H. armigera during an early stage of HearNPV infection. Future studies including RNAi, over expression or protein manipulation are still needed for a deeper understanding of the identified DEGs and will shed light on the interaction between H. armigera and HearNPV.

Acknowledgements This study was supported by grants from the National ‘‘973’’ Project (No. 2012CB114103). The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Appendix A. Supplementary material Fig S1. Sequencing saturation analysis of two DGEs.

19

Fig S2. Differentially expressed genes in control vs HearNPV-inoculated libraries. Table S1. Information of primers used for qRT-PCR. Table S2. Percentage of clean reads with RPKM value in the specific range.

References Audic S., Claverie J.M., 1997. The significance of digital gene expression profiles. Genome Res. 7, 986−995. Bao Y., Tang X., Lv Z., Wang X., Tian C., Xu Y., Zhang C., 2009. Gene expression profiling of resistant and susceptible Bombyx mori strains reveals nucleopolyhedrovirus-associated variations in host gene transcript levels. Genomics. 94, 138−145. Bao Y., Lv Z., Liu Z., Xue J., Xu Y., Zhang C., 2010. Comparative analysis of Bombyx mori nucleopolyhedrovirus responsive genes in fat body and haemocyte of B. mori resistant and susceptible strains. Insect Mol. Biol. 19, 347−358. Bennett S.T., Barnes C., Cox A., Davies L., Brown C., 2005. Toward the $1000 human genome. Pharmacogenomics. 6, 373−382. Breitenbach J.E., Popham J.R., 2013. Baculovirus replication induces the expression of heat shock proteins in vivo and in vitro. Arch. Virol. 158, 1517−1522. Chen X., Ijkel W.F.J., Tarchini R., Sun X., Sandbrink H., Wang H., Peters S., Zuidema D., Lankhorst R.K., Vlak J.M., Hu Z., 2001. The sequence of the Helicoverpa armigera single nucleocapsid nucleopolyhedrovirus genome. J. Gen. Virol. 82, 241–257. Chen Y., Zhong S., Fei Z., Gao S., Zhang S., Li Z., Wang P., Blissard G.W., 2014. Transcriptome response of the host Trichoplusia ni to infection by the baculovirus Autographa californica multiple nucleopolydedrovirus. J. Viol. 88, 13781–13797. Choi J.Y., Roh J.Y., Wang Y., Tao X.Y., Lee J.H., Liu Q., Kim J.S., Shin S.W., Je Y.H., 2012. Analysis of genes expression of Spodoptera exigua larvae upon AcMNPV infection. PLoS ONE. 7, e42462. Etebari K., Matindoost L., Mirhoseini S.Z., Turnbull M.W., 2007. The effect of BmNPV infection on protein metabolism in silkworm (Bombyx mori) larva. Invert. Surviv. J. 4, 13−17. Fan X., Lu M., Meng X., Rui C., 2003. An improved rearing technique for Helicoverpa armigera. 20

Chinese Bull. Entomol. 40, 85−87. (in Chinese) Gao K., Deng X., Qian H., Qin G., Guo X., 2014. Digital gene expression analysis in the midgut of 4008 silkworm strain infected with cytoplasmic polyhedrosis virus. J. Invertebr. Pathol. 115, 8−13. Gatehouse H.S., Poulton J., Markwick N.P., Gatehouse L.N., Ward V.K., Yong V.L., Luo Z., Schaffer R., Christeller J.T., 2009. Changes in gene expression in the permissive larval host lightbrown apple moth (Epiphyas postvittana, Tortricidae) in response to EppoNPV (Baculoviridae) infection. Insect Mol. Biol. 18, 635–648. Georgievska L., Joosten N., Hoover K., Cory J.S., Vlak J.M., van der Werf W., 2010. Effects of single and mixed infections with wild type and genetically modified Helicoverpa armigera nucleopolyhedrovirus on movement behaviour of cotton bollworm larvae. Entomol. Exp. Appl. 135, 56–67. Grzywacz D., Rabindra R.J., Brown M., Jones K.A., Parnell M., 2000. The Helicoverpa armigera NPV production manual. http://www.fao.org/docs/eims/upload/agrotech. Hawtin R.E., Zarkowska T., Amold K., Thomas C.J., Gooday G.W., King L.A., Kuzio J.A., Possee R.D., 1997. Liquefaction of Autographa californica multiple nuclear polyhodrosis virus infected insects is dependent on the integrity of virus encoded chitinase and cathepsin genes. Virology. 238, 243−253. Hoover K., Grove M., Gardner M., Hughes D.P., McNeil J., Slavicek J., 2011. A gene for an extended phenotype. Science. 333, 1401. Hou C., Qin G., Liu T., Geng T., Gao K., Pan Z., Qian H., Guo X., 2014. Transcriptome analysis of silkworm, Bombyx mori, during early response to Beauveria bassiana challenges. PLoS ONE. 9, e91189. Iatrou K., Ito K., Witkiewicz H., 1985. Polyhedrin gene of Bombyx mori nuclear polyhedrosis virus. J. Virol. 54, 436−445. Iwanaga M., Shibano Y., Ohsawa T., Fujita T., Katsuma S., Kawasaki H., 2014. Involvement of HSC70-4 and other inducible HSPs in Bombyx mori nucleopolyhedrovirus infection. Virus Res. 179, 113−118. Kamita S.G., Nagasaka K., Chua J.W., Shimada T., Mita K., Kobayashi M., Maeda S., Hammock 21

B.D., 2005. A baculovirus-encoded protein tyrosine phosphatase gene induces enhanced locomotory activity in a lepidopteran host. Proc. Natl. Acad. Sci. U S A. 102, 2584−2589. Kanehisa M., Araki M., Goto S., Hattori M., Hirakawa M., Itoh M., Katayama T., Kawashima S., Okuda S., Tokimatsu T., Yamanishi Y., 2008. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, D480−484. Lanier L.M., Slack J.M., Volkman L.E., 1996. Actin binding and proteolysis by the baculovirus AcMNPV: the role of virion-associated V-CATH. Virology. 216, 380−388. Li R., Yu C., Li Y., Lam T.W., Yiu S.M., Kristiansen K., Wang J., 2009. SOAP2: An improved ultrafast tool for short read alignment. Bioinformatics. 25, 1966−1967. Maynard N.D., Gutschow M.V., Birch E.W., Covert M.W., 2010. The virus as metabolic engineer. Biotechnol. J. 5, 686−694. Monastirioti M., 1999. Biogenic amine systems in the fruit fly Drosophila melanogaster. Microsc. Res. Tech. 45, 106−121. Monteiro F., Carinhas N., Carrondo M.J., Bernal V., Alves P.M., 2012. Toward system-level understanding of baculovirus-host cell interactions: from molecular fundamental studies to large-scale proteomics approaches. Frontier Microbiol. 3, 9. Morozova O., Marra M.A., 2008. Applications of next-generation sequencing technologies in functional genomics. Genomics. 92, 255−264. Mortazavi A., Williams B.A., McCue K., Schaeffer L., World B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods. 5, 621−628. Moscardi F., 1999. Assessment of the application of Baculoviruses for control of Lepidoptera. Annu. Rev. Entomol. 44, 257–289. Nguyen Q., Palfreyman R.W., Chan L.C.L., Reid S., Nielsen L.K., 2012. Transcriptome sequencing of and microarray development for a Helicoverpa zea cell line to investigate in vitro insect cell-baculovirus interactions. PLoS ONE. 7, e36324. Nguyen Q., Chan L.C.L., Nielsen L.K., Reid S., 2013. Genome scale analysis of differential mRNA expression of Helicoverpa zea insect cells infected with a H. armigera baculovirus. Virology. 444, 158−170. Nobiron I., O’Reilly D.R., Olszewski J.A., 2003. Autographa californica nucleopolyhedrovirus 22

infection of Spodoptera frugiperda cells: a global analysis of host gene regulation during infection, using a differential display approach. J. Gen. Virol. 84, 3029–3039. Okano K., Shimada T., Mita K., Maeda S., 2001. Comparative expressed-sequence tag analysis of differential gene expression profiles in BmNPV-infected BmN cells. Virology. 282, 348–356. Ooi B.G., Miller L.K., 1988. Regulation of host RNA levels during baculovirus infection. Virology. 166, 515−523. O’Reilly D.R., Brown M.R., Miller L.K., 1992. Alteration of ecdysteroid metabolism due to baculovirus infection of the fall armyworm Spodoptera frugiperda: Host ecdysteroids are conjugated with galactose. Insect Biochem. Mol. Biol. 22, 313−320. Popham H.J.R., Grasela J.J., Goodman C.L., McIntosh A.H., 2010. Baculovirus infection influences host protein expression in two established insect cell lines. J. insect physiol. 56, 1237−1245. Riegel C.I., Slavicek J.M., 1997. Characterization of the replication cycle of the Lymantria dispar nuclear polyhedrosis virus. Virus Res. 51, 9−17. Sagisaka A., Fujita K., Nakamura Y., Ishibashi J., Noda H., Imanishi S., Mita K., Yamakawa M., Tanaka H., 2010. Genome-wide analysis of host gene expression in the silkworm cells infected with Bomyx mori nucleopolyhedrovirus. Virus Res. 147, 166–175. Salem T.Z., Zhang F., Xie Y., Thiem S.M., 2011. Comprehensive analysis of host gene expression in Autographa californica nucleopolyhedrovirus-infected Spodoptera frugiperda cells. Virology. 412, 167–178. Thomson E., Ferreira-Cerca S., Hurt E., 2013. Eukaryotic ribosome biogenesis at a glance. J. Cell Sci. 126, 4815−4821. Yang X., Liu X., Xu X., Li Z., Li Y., Song D., Yu T., Zhu F., Zhang Q., Zhou X., 2014. Gene expression profiling in winged and wingless cotton aphids, Aphis gossypii (Hemiptera: Aphididae). Int. J. Biol. Sci. 10, 257−267. Yue Y., Tang X., Xu L., Yan W., Li Q., Xiao S., Fu X., Wang W., Li N., Shen Z., 2015. Early responses of silkworm midgut to microsporidium infection-digital gene expression analysis. J. Invertebr. Pathol. 124, 6−14.

23

Figure Captions Fig.1 Transcriptional profile of NPV-encoded polh in different tissues of H. armigera larvae post inoculation with 105 OBs of HearNPV. Relative gene expressional level of polh in the head (H), midgut (MG), ventral nerve cord (VNG), fat body (FB) and Malpigian tube (MT) were tested at 24, 48, 72 and 96 hpi using qRT-PCR. Fig.2 Comparison analysis of qRT-PCR detection with corresponding DEG profiling. The y-axis indicates the relative transcriptional level of genes quantified by fold change (log2Ratio) between HearNPV-inoculated and control samples at 24 hpi. The x-axis represents the specific assessed genes, as 1: carboxyl/choline esterase, 2: aldo-keto reductase, 3: fatty acid-binding protein 1, 4: alpha amylase, 5: ecdysteroid-regulated protein, 6: neutral lipase, 7: astacin, 8: glucose oxidase-like enzyme, 9: fructosidase, 10: pancreatic lipase 2, 11: trypsin, 12: juvenile hormone diol kinase, 13: uridine 5'-monophosphate synthase, 14: stress-induced phosphoprotein 1, 15: heat shock protein 70, 16: serine protease 52, 17: heat shock protein 90, 18: ribonucleoside-diphosphate reductase, 19: carboxypeptidase precursor, 20: palmitoyl-protein thioesterase 1. Fig.3 Time course validation about transcriptional level of DEGs using qRT-PCR. Genes related to A, defense response (stress-induced phosphoprotein 1, heat shock protein 70, and serine protease 52); B, salivary secreted enzymes (glucose oxidase-like enzyme, fructosidase, carboxyl/choline esterase); C, nutrition metabolism (trypsin and carboxypeptidase precursor, and neutral lipase); and D, regulation of growth and development (ecdysteroid-regulated protein and juvenile hormone diol kinase) were detected at 0, 6, 12, 24, 48, 72, 96 hpi. The y-axis indicates the relative transcriptional level of tested genes related to gene expression level at 0 hpi. The x-axis represents the tested time points post HearNPV inoculation. “*” on bar represents the significant difference of gene expression level between HearNPV-inoculated and control larvae at P<0.05, and “**” represents the extreme significance at P<0.01.

Fig.4 DEG network. A, potential gene network analysis of the up-regulated DEGs (∣log2 Ratio∣≥2 ), the genes in circles of dark red represent the up-regulated genes identified in present DGE analysis, the genes in pink circles represent other related genes in the network predicted by GeneMANIA; B, potential gene network analysis of the down-regulated DEGs (∣log2 Ratio∣≥2 ), the genes in circles of dark green represent the down-regulated genes identified in present DGE analysis, the genes in circles of light green represent other related genes in the network predicted by GeneMANIA.

24

47500

Relative expression level of polh

45000 42500 40000

H MG VNG FB MT

a

37500 35000

b

3000 2500

a

2000 1500

b

1000 500

c d a bc

0 24

48

c c

bc

bc c

c 72

Post NPV inoclulation (hrs)

96

25 Gene expression in DGE Gene expression in qRT-PCR

Relative expression level (Fold change=log 2 Ratio)

20 15 10 5 0 -5 -10 -15 -20 1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20

Gene number

carboxyl/choline esterase stress-induced phosphoprotein 1

NPV

40.0

16

**

5

**

75

**

4

14

70

3

** **

10 8 6

1

**

4

**

2 0 0

0

B

**

6

6

1212

2424

4848

*

**

6

12

24

48

30

16 6

20

** **

2 0 0

6

12

24

48

72

96

ddH2O HaSNPV

12

**

**

24

48

72

96

carboxyl/choline esterase

** **

75

** 70

** ** * 0

96

0

80

fructosidase

12 10 8 6 4 2 0

**

**

72

40

17

4

6

**

.1

50

18

**

0.0 0

60

19

**

2

96 96

glucose oxidase-like enzyme

1.0

**

0

7272

10.0 .2

4

**

**

20.0

10 6

** **

serine protease 52

30.0

**

12

2

**

**

12

0

Relative gene expression level

heat shock protein 70

ddH2O

80

Relative gene expression level

Relative gene expression level

A

6

12

24

48

72

**

12 10 8 6 4 2 0

**

** ** 0

96

6

12

24

48

72

96

Relative gene expression level

C 30

trypsin

**

8

carboxypeptidase precursor

7

**

20

**

6

10

**

2

**

2

**

2

1

6

12

24

48

72

96

**

**

12

24

48

0

0

0

** **

1

**

0

**

3

** *

6

4

**

4

**

**

5

**

3

neutral lipase

0

6

12

24

48

72

96

0

6

72

96

Time post NPV inoculation (hours)

Relative gene expression level

D ecdysteroid-regulated protein

juvenile hormone diol kinase 28

60

**

**

55

26

50

**

25

12 10 8 6 4 2 0

**

20 15

**

10 5

**

*

0 0

6

12

24

48

72

96

Time post NPV inoculation (hours)

** **

0

6

12

24

**

48

72

96

Time post NPV inoculation (hours)

A

B

Table 1. Distributions of sequencing reads in two DGE libraries Category

Parameter

Reference gene

Control

HearNPV-inoculated

ESTs of H. armigera

Genome of B. mori

ESTs of H. armigera

Genome of B. mori

Raw data

total reads

4842523

4842523

4765688

4765688

Clean reads

total reads

4817921

4817921

4736444

4736444

total basepair

236098129

236098129

232085756

232085756

total mapped reads

1261224

143836

2050163

253706

percentage to reference(%)

26.18

2.99

43.28

5.36

total mapped reads

1074667

121749

1763991

200548

percentage to reference(%)

22.31

2.53

37.24

4.47

total unmapped reads

3556697

4674085

2686281

4482738

percentage to reference(%)

73.82

97.01

56.72

94.64

All reads mapping to gene

Unique reads mapping to gene

Unmapped reads

25

Table 2. Gene categories of DEGs. Up-regulated Gene categories

Down-regulated

No. of gene

Percent %

No. of gene

Percent %

47

34.56

191

31.99

metabolic process

34

25

56

9.38

cellular process

2

1.47

60

10.05

biological regulation

5

3.68

46

7.71

translation

0

0

14

2.35

defense response

6

4.41

15

2.51

Biological process

Molecular function

68

50

237

39.7

catalytic activity

39

28.68

40

6.7

binding

12

8.82

157

26.3

enzyme regulator activity

13

9.56

37

6.2

transporter activity

4

2.94

3

0.5

21

15.44

169

28.31

Cellular component

26

Table3. Significantly up-regulated genes in H. armigera larvae at 24h post-inoculation of HearNPV. Gene ID

Fold change§

FDR*

Homologs (Acession No.)

Score

E-value

354972604

18.74

0

Helicoverpa armigera, carboxyl/choline esterase (FJ997317)

1358

0

354968175

14.15

9.16E-05

Bombyx mori, indole-3-acetaldehyde oxidase-like (XM_004928315)

100

2.00E-17

354969115

4.22

0.00042617

Helicoverpa armigera, D-type cecropins like (EU041763)

482

2.00E-132

354969741

4.07

6.58E-05

Agrotis ipsilon, aldo-keto reductase (KC007347)

536

1.00E-148

354980173

4.02

4.46E-12

Microplitis demolitor, probable phytanoyl-CoA dioxygenase (XM_008550939)

425

4.00E-115

194240282

3.98

2.17E-07

Helicoverpa armigera, fatty acid-binding protein 1 (EU325558)

120

5.00E-24

112350170

3.90

1.17E-05

Bombyx mori, alpha amylase (NM_001195462 )

237

2.00E-58

112349948

3.82

5.53E-12

Helicoverpa armigera, ecdysteroid-regulated protein (DQ875272)

544

7.00E-151

354979253

3.60

0

Helicoverpa armigera, cuticle protein 1 (EU526836)

778

0.00E+00

112349806

3.56

1.33E-09

Helicoverpa armigera, neutral lipase (JF999960)

1361

0.00E+00

112350220

3.54

1.42E-12

Helicoverpa armigera, peptidoglycan recognition protein B (PGRP-B) (JX082166)

917

0.00E+00

194240091

3.48

2.04E-13

Bombyx mori, cuticular protein RR-1 motif 41 (NM_001044025)

91.5

3.00E-15

354970709

3.14

1.59E-10

Helicoverpa armigera, glucose oxidase-like enzyme (EU629216)

1296

0

354972955

3.38

2.16E-07

Helicoverpa armigera, attacin (AY948540)

1021

0

354965775

3.37

0.0001

Helicoverpa armigera, cytochrome P450 CYP4L5 (KM016727)

432

2.00E-117

112350232

3.32

7.30E-09

Spodoptera frugiperda, astacin (EU915282 )

509

1.00E-140

354968676

3.14

0.00098154

Bombyx mori, putative inorganic phosphate cotransporter-like (XM_004925796)

455

2.00E-124

354966822

2.73

1.17E-05

Helicoverpa armigera, fructosidase (EF600050)

1240

0.00E+00

354971081

2.69

1.17E-05

Helicoverpa armigera, lysozyme (HM588761)

776

0.00E+00

112349856

2.69

7.29E-09

Mamestra configurata, pancreatic lipase 2 (EU660854)

670

0.00E+00

194240144

2.66

4.95E-07

Helicoverpa armigera, chitin deacetylase 5b (cda5b) (GQ411191)

289

5.00E-75

354963324

2.65

1.18E-05

Helicoverpa armigera, gloverin (Glo) (GU182908)

623

7.00E-175

27

112349802

2.59

6.60E-12

Helicoverpa armigera, trypsin (EU325548)

1366

0.00E+00

354965313

2.47

1.12E-06

Bombyx mori, sepiapterin reductase Spr (AB465551 )

114

1.00E-21

354970948

2.38

1.73E-12

Helicoverpa zea, lysozyme precursor (FJ535250)

1142

0.00E+00

112349782

2.36

2.03E-12

Plutella xylostella, Juvenile hormone diol kinase (AB180433 )

251

6.00E-63

354963507

2.18

1.82E-12

Bombyx mori, pyrroline-5-carboxylate reductase-like (XM_004926811)

242

2.00E-60

354962641

2.06

2.55E-12

Bombyx mori, hydroxyacid oxidase 1-like (XM_004924614)

370

8.00E-99

Note:

§

Fold change represents log2 Ratio; Ratio=expression level of genes in HearNPV infected larvae/ expression level of genes in control larvae, at 24h post HearNPV inoculation; expression level of genes were calculated with algorithm of RPKM (Reads Per Kb per Million reads) (Mortazavi et al., 2008).

*

FDR (False Discovery Rate) ≤0.001 was used as the threshold to judge the significance of gene expression difference (Benjamini et al., 2001).

28

Table 4. Significantly down-regulated genes in H. armigera larvae at 24h post-inoculation of HearNPV. Gene ID

Fold change§

FDR*

Homologs( Acession No.)

Score

E-value

354978081

-14.98

3.92E-10

Papilio xuthus, palmitoyl-protein thioesterase 1 (AK401487)

190

2.00E-44

354962003

-13.83

0.00031

Bombyx mori, probable dolichyl pyrophosphate Glc1Man9GlcNAc2

156

3.00E-34

237

1.00E-58

alpha-1,3-glucosyltransferase-like (XM_004927809) 354962136

-13.55

0.00074

Bombyx mori, UDP-xylose and UDP-N-acetylglucosamine transporter-like (XM_004928594)

354977630

-13.42

0.00074

Saccoglossus kowalevskii, glucose dehydrogenase [FAD, quinone]-like (XM_002731972)

73.4

4.00E-09

112350216

-8.19

7.99E-72

Spodoptera frugiperda, ribosomal protein S18（AF400215）

650

0

354980754

-6.47

2.15E-20

Antheraea mylitta, retrotransposon polyprotein gene （AF530470）

230

2.00E-56

354979323

-5.97

5.10E-71

Helicoverpa armigera, small heat shock protein G16 (KC689795)

1321

0

354980773

-5.28

1.46E-91

Helicoverpa armigera, small heat shock protein G8 (KC689794)

1027

0

354969241

-4.62

3.27E-05

Bombyx mori,condensin complex subunit 2-like （XM_004931542）

369

3.00E-98

354968334

-4.30

0.00039

Bombyx mori,nudC domain-containing protein 1-like （XM_004928650）

176

3.00E-40

354981978

-4.30

0.00039

Apis florea, c-Myc-binding protein-like (XM_003689737)

111

9.00E-21

354964458

-4.17

5.92E-07

Oryzias latipes,nardilysin-like （XM_004082220）

91.5

1.00E-14

354978025

-3.88

1.51E-05

Helicoverpa armigera,cytochrome P450 (CYP337B3v1) (JQ995292)

930

0.00E+00

354968311

-3.58

2.82E-27

Spodoptera frugiperda, DnaJ (KF562156)

720

0

354963107

-3.35

4.01E-08

Bombyx mori, leucine-rich repeat-containing protein 47-like (XM_004928928)

260

1.00E-65

354973166

-3.33

2.25E-12

Bombyx mori, serine/threonine-protein kinase RIO2-like （XM_004925387）

297

1.00E-76

354970652

-3.24

5.11E-06

Bombyx mori, nucleosome-remodeling factor subunit NURF301-like (XM_004924801)

340

1.00E-89

354968380

-3.21

4.30E-11

Bombyx mori,nucleolar protein 14 homolog （XM_004932959）

535

3.00E-148

22474353

-3.17

2.98E-35

Bombyx mori, actin cytoskeleton-regulatory complex protein PAN1-like

95.1

9.00E-16

（XM_004922373） 354962622

-2.91

2.39E-07

Bombyx mori, large neutral amino acids transporter small subunit 1-like

489

1.00E-134

(XM_004930918)

29

354970590

-2.79

7.03E-05

Bombyx mori, nitric oxide synthase-interacting protein （XM_004927696）

298

3.00E-77

22474248

-2.76

1.65E-238

Helicoverpa armigera, carboxypeptidase precursor (AJ626862)

1356

0.00E+00

354973308

-2.76

3.52E-07

Bombyx mori, ATPase WRNIP1-like （XR_209986）

135

8.00E-28

354966890

-2.67

1.44E-06

Papilio polytes, septin interacting protein 1 (AK405351)

365

4.00E-97

354980520

-2.65

6.89E-13

Bombyx mori, DNA polymerase delta subunit 2-like (XM_004931391)

361

3.00E-96

354978843

-2.63

6.74E-25

Bombyx mori, ribonucleoside-diphosphate reductase large subunit-like (XM_004925495)

607

7.00E-170

354970365

-2.52

0.00019

Bombyx mori, midasin-like (XM_004922432)

298

5.00E-77

354978803

-2.51

1.21E-08

Bombyx mori,leucine--tRNA ligase （XM_004927710）

351

6.00E-93

354982429

-2.48

1.30E-06

Bombyx mori, RNA-binding protein 5-like (XR_209675)

258

7.00E-65

354966244

-2.44

2.51E-06

Bombyx mori, cell cycle control protein 50A-like (XM_004928819)

318

4.00E-83

354974310

-2.42

2.85E-07

Bombyx mori, proteasome activator complex subunit 3-like (XM_004932201)

320

1.00E-83

354965623

-2.42

4.18E-05

Bombyx mori, protein Peter pan-like (XM_004925204)

260

2.00E-65

354962651

-2.40

1.43E-05

Bombyx mori, isoleucine--tRNA ligase, cytoplasmic-like (XM_004929636)

172

4.00E-39

354978874

-2.36

8.14E-05

Bombyx mori, MYST histone acetyltransferase (mof) (DQ442997)

590

5.00E-165

354980786

-2.35

9.29E-06

Bombyx mori, talin-2-like (XR_209778)

360

1.00E-95

354974921

-2.30

3.57E-151

Helicoverpa armigera, heat shock protein (JF417984)

1038

0.00E+00

112349947

-2.30

6.97E-21

Mamestra configurata, serine protease 52 (HM990174)

204

1.00E-48

22474229

-2.26

5.25E-44

Bombyx mori, ribonucleoside diphosphate reductase small subunit (RnrS)

400

1.00E-107

(NM_001257337) 354974061

-2.26

3.85E-06

Bombyx mori, choline/ethanolaminephosphotransferase 1-like (XM_004927681)

392

2.00E-105

354975382

-2.25

1.88E-30

Erinaceus europaeus,stress-induced phosphoprotein 1 (STIP1) （XM_007518548）

73.7

3.00E-09

22474295

-2.24

0.00029

Agrotis segetum, putative peroxisomal acyl-CoA oxidase 1(KJ622089)

331

5.00E-87

354973563

-2.23

0.00087

Bombyx mori,protein Mpv17-like (XM_004928465)

91.1

1.00E-14

354973669

-2.20

5.34E-07

Apis dorsata, serine/threonine-protein phosphatase 1 regulatory subunit 10-like

77.6

2.00E-10

266

2.00E-67

(XM_006613828) 354966230

-2.17

0.00055

Bombyx mori, adaptor protein AP-2 complex subunit alpha-like (XM_004932091)

30

354963680

-2.17

1.35E-05

Bombyx mori, probable ATP-dependent RNA helicase Dbp45A-like

283

2.00E-72

(XM_004922618) 354979834

-2.17

1.00E-06

Bombyx mori, uridine 5'-monophosphate synthase (NM_001046695)

162

6.00E-36

354981183

-2.16

3.97E-05

Bombyx mori, oligoribonuclease, mitochondrial-like (XM_004932658）

69.9

3.00E-08

354979353

-2.15

3.26E-10

Bombyx mori, putative ribosomal RNA methyltransferase NOP2-like (XM_004930390)

185

8.00E-43

354973784

-2.13

2.47E-05

Bombyx mori, DNA supercoiling factor (JN171664)

248

4.00E-62

194240237

-2.11

1.17E-09

Bombyx mori, serine protease inhibitor 10 (serpin-10)(NM_001146231)

244

6.00E-61

354967621

-2.09

1.52E-10

Bombyx mori, deoxynucleoside kinase-like (NM_001043554)

282

5.00E-72

354974438

-2.08

9.89E-06

Bombyx mori, chromatin-remodeling complex ATPase chain Iswi-like (XM_004927827)

467

6.00E-128

354974515

-2.07

3.54E-28

Bombyx mori, mitochondrial import inner membrane translocase (NM_001098371)

210

2.00E-50

354968642

-2.07

0.00063

Bombyx mori, ESF1 homolog (XM_004932024)

269

2.00E-68

354980760

-2.05

1.34E-06

Bombyx mori, elongation factor G, mitochondrial-like (XM_004929292)

285

5.00E-73

354970588

-2.04

8.36E-07

Bombyx mori, large subunit GTPase 1 homolog (XM_004932875)

192

3.00E-45

354973056

-2.04

0.00040

Spodoptera exigua, topoisomerase 1B (JN258956)

485

3.00E-133

354980564

-2.04

3.46E-12

Bombyx mori, ubiquitin-conjugating enzyme E2 C-like (XM_004925988)

497

9.00E-137

354963105

-2.01

2.86E-09

Bombyx mori,splicing factor 3B subunit 1-like (XM_004924775)

627

6.00E-176

354966977

-2.01

2.91E-22

Heliconius melpomene, DEAD box ATP-dependent RNA helicase-like protein (GQ452007)

623

9.00E-175

354978338

-2.01

6.65E-08

Bombyx mori, microfibrillar-associated protein 1-like （XM_004929597）

179

3.00E-41

354962771

-2.00

2.46E-12

Bombyx mori, transcription factor Dp-1-like (XM_004927623)

433

1.00E-117

Note:

§

Fold change represents log2 Ratio; Ratio=expression level of genes in HearNPV infected larvae/ expression level of genes in control larvae, at 24h post HearNPV inoculation; expression level of genes were calculated with algorithm of RPKM (Reads Per Kb per Million reads) (Mortazavi et al., 2008).

*

FDR (False Discovery Rate) ≤0.001 was used as the threshold to judge the significance of gene expression difference (Benjamini et al., 2001).

31

Table 5. Significant KEGG pathways at 24h post-inoculation of HearNPV KEGG pathways

Pathway ID

P value*

No. of gene up-regulated

down-regulated

Ribosome biogenesis in eukaryotes

ko03008

0.000256

31

Prion diseases

ko05020

0.000637

9

Pancreatic secretion

ko04972

0.001118

Aminoacyl-tRNA biosynthesis

ko00970

0.004548

Protein processing in endoplasmic reticulum

ko04141

0.004839

8

23

Protein digestion and absorption

ko04974

0.011471

14

6

Pyrimidine metabolism

ko00240

0.013527

1

17

Salivary secretion

ko04970

0.020611

5

1

Valine, leucine and isoleucine biosynthesis

ko00290

0.027449

6

Spliceosome

ko03040

0.02895

24

p53 signaling pathway

ko04115

0.036668

Antigen processing and presentation

ko04612

0.038324

Neuroactive ligand-receptor interaction

ko04080

0.049727

21

6 19

1

3 10

5

7

*

Note: Pathways with P value ≤0.05 are regarded significantly enriched in DEGs (Differentially expressed genes).

32

Graphical Abstract

A

C

47500

Relative expression level of polh

45000 42500 40000

up-regulated

H MG VNG FB MT

25% metabolic process

2.94%

cellular process

37500 35000

biological regulation

9.56%

b

1.47%

translation

3.68%

defense response

3000 2500

a

2000

catalytic activity

1500

b

1000 500

c d a bc

0 24

48

c c

bc

4.41%

transporter activity 72

96

down-regulated

9.38%

28.31%

up regulated down regulated

cellular component

28.68%

10.05%

metabolic process

7.71%

Number of DEGs

binding enzyme regulator activity

c

c

700

600

8.82%

bc

Post NPV inoclulation (hrs)

B

15.44%

a

500

cellular process biological regulation tanslation

400

defense response 300

catalytic activity binding

200

enzyme regulator activity transporter activity

100

2.35% 6.70% 2.51%

6.20%

0.50%

0

Control vs NPV inoculated for 24h

26.30%

cellular component

Highlights  Transcriptional profile was analyzed at 24 h post HearNPV infection via DGE.  A total of 136 up-regulated and 597 down-regulated DEGs were identified.  Random selected DEGs were revalidated by qRT-PCR.  Increased metabolism and hormone disorder occurred in HearNPV-infected larvae.  The results provide basic data for the molecular mechanism of HearNPV infection.

33