key: cord-0002218-z1b7pj3s
authors: Wang, Xue-Yang; Yu, Hai-Zhong; Geng, Lei; Xu, Jia-Ping; Yu, Dong; Zhang, Shang-Zhi; Ma, Yan; Fei, Dong-Qiong
title: Comparative Transcriptome Analysis of Bombyx mori (Lepidoptera) Larval Midgut Response to BmNPV in Susceptible and Near-Isogenic Resistant Strains
date: 2016-05-11
journal: PLoS One
DOI: 10.1371/journal.pone.0155341
sha: 345a04a8633d593b5a400b559d533a6b4595f490
doc_id: 2218
cord_uid: z1b7pj3s

Bombyx mori nucleopolyhedrovirus (BmNPV) is one of the primary pathogens causing severe economic losses in sericulture. However, the molecular mechanism of silkworm resistance to BmNPV remains largely unknown. Here, the recurrent parent P50 (susceptible strain) and the near-isogenic line BC9 (resistance strain) were used in a comparative transcriptome study examining the response to infection with BmNPV. A total of 14,300 unigenes were obtained from two different resistant strains; of these, 869 differentially expressed genes (DEGs) were identified after comparing the four transcriptomes. Many DEGs associated with protein metabolism, cytoskeleton, and apoptosis may be involved in the host response to BmNPV infection. Moreover, some immunity related genes were also altered following BmNPV infection. Specifically, after removing genetic background and individual immune stress response genes, 22 genes were found to be potentially involved in repressing BmNPV infection. These genes were related to transport, virus replication, intracellular innate immune, and apoptosis. Our study provided an overview of the molecular mechanism of silkworm resistance to BmNPV infection and laid a foundation for controlling BmNPV in the future.

The silkworm, Bombyx mori L. (Lepidoptera: Bombycidae) has been domesticated for production of cocoons for more than 5000 years. Silkworm rearing and the silk industry still play an important role in China, India and many other developing countries. B. mori is also a good model for the study of insect genetics and immunology [1] [2] [3] [4] . Bombyx mori nucleopolyhedrovirus (BmNPV) is the principal silkworm pathogen and causes serious economic losses in sericulture every year. Among numerous silkworm strains, most are susceptible to BmNPV infection, although a few strains exhibit resistance [5] . The heredity of silkworm resistance against BmNPV infection is a relatively complicated process because resistance is controlled both by major dominant genes and multiple genes of micro-effect [6] .

A series of studies have made significant progress in understanding silkworm resistance against BmNPV infection. Xu et al. reported that Bms3a was potentially involved in resistance to BmNPV infection [7, 8] . B. mori lipase-1, serine protease-2 and alkaline trypsin protein extracted from the digestive juice of larvae midguts showed strong antiviral activity in vitro [9] [10] [11] . Using comparative proteomics, arginine kinase was found to be involved in the antiviral process of different resistant strains of silkworm [12] . In our laboratory, a total of 12 proteins that are potentially involved in viral infection were identified using one-and two-dimensional electrophoresis followed by virus overlay assays. These proteins could be categorized into the following groups: endocytosis, intracellular transportation, and host responses [13] . Immune responses were found to be synergistically regulated by the Toll, Janus kinase/signal transducers and activators of the transcription (JAK/STAT) and immune deficiency (IMD) pathways, which could act as an important defense against exogenous pathogenic infection in conjunction with subsistent pathogen recognition receptors and response proteins [14] [15] [16] [17] [18] . However, the molecular mechanisms of silkworm resistance to BmNPV infection are still not fully elucidated.

In recent years, the high-throughput nature of next generation sequencing (NGS), using platforms such as Illumina HiSeq TM 2500 have provided fascinating opportunities in the life sciences and dramatically improved the efficiency and speed of gene discovery, especially in the research of host cell responses to exogenous pathogenic infection [19] . For example, Hu et al. obtained numerous differentially expressed genes (DEGs) involved in metabolism, immunity, and inflammatory responses in Microtus fortis following infection with Schistosoma japonicum based on comparative transcriptome analysis [20] . Diege et al. examined different fish tissues infected with salmon anemia virus (ISAV) using high-throughput transcriptomics and found a strong correlation between functional modules and viral-segment transcription [16] . NGS technology was also used to explore the molecular mechanism of silkworm resistance against exogenous pathogens. Kolliopoulou et al. reported that several genes related to physical barriers, immune response, proteolytic/metabolic enzymes, heat-shock proteins, and hormonal signaling were possibly involved in silkworm resistance against Bombyx mori cytoplasmic polyhedrosis virus (BmCPV) infection; although these genes might be induced by the virus in order to increase infectivity [21] . Additionally, several candidate genes, such as BmEts, BmToll10-3 and Hsp20-1, have been identified in the initial stage of BmNPV infection by analyzing the global transcriptional profile of silkworm cell lines and heads following BmNPV infection [22, 23] .

In order to gain a global view of the molecular changes in silkworms during BmNPV infection, we selected near-isogenic line BC9 and recurrent parent P50 for transcriptome sequencing. Through comparative analysis of the transcriptomes from these two strains, a total of 869 DEGs were obtained, which included many genes potentially related to BmNPV-resistance. Our results may provide some reliable evidence to clarify the BmNPV-resistance molecular mechanism in silkworm.

BmNPV (T3 strain) was maintained in the Key Laboratory of Sericulture, Anhui Agricultural University, Hefei, China. The virus was obtained from the hemolymph of infected larvae and purified by repeated and differential centrifugation according to the protocol developed by Rahman et al. [24] . The concentration of the virus (OB/mL) was determined by hemocytometer.

The recurrent parent P50 (susceptible strain), the donor parent A35 (resistant strain) and the near-isogenic line BC9 were maintained in our laboratory. The near-isogenic line was constructed according to the protocol used by Yao et al. [6] . In brief, the recurrent parent P50 was crossed to the donor parent A35; progeny were repeatedly backcrossed with the recurrent parent for nine generations and each progeny was screened by BmNPV.

The first three instar larvae were reared on a fresh artificial diet at 26±1°C, 75±5% relative humidity, and a 12 hours day/night cycle. The rearing temperature for the last two instars was reduced to 24±1°C; other conditions were unchanged. On the first day of fifth instar, all larvae were starved for 24 hours and then fed with 5 μL BmNPV suspended in sterile water (1.0×10 5 OB/mL) per larva orally; the control group was treated with sterile water. BmNPV occlusion bodies (OB) began fast proliferation at approximately 24 hours post inoculation (hpi) [25] ; thus, this time was considered optimal for sample collection. Silkworm larvae were dissected and the midgut tissues were removed and then washed in PBS (137 mM NaCl, 2.7 mM KCl, 4.3 mM Na 2 HPO 4 , and 1.4 mM KH 2 PO 4 , pH 7.4) prepared with diethy pyrocarbonate (DEPC) (Sangon, China) treated H 2 O. Thirty larvae midguts were mixed together to minimize individual genetic differences. Samples were flash-frozen in liquid nitrogen and pulverized, and 100 mg of sample were added directly into a RNAase free microcentrifuge tube containing 1.0 mL TRIzol Reagent (Invitrogen, USA) and stored at -80°C for later use.

The level of silkworm resistance to BmNPV was tested following the protocol developed by Cheng et al. [26] . The fourth instar larvae were inoculated with BmNPV at different concentrations; inoculations were conducted in triplicate. The level of silkworm resistance was calculated using IBM SPSS Statistics 20 (IBM, USA).

The silkworm midguts dissolved in TRIzol Reagent were homogenized. Total RNA of midguts were extracted according to the manufacturer's protocol. Concentrations were quantified using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). The purity of all RNA samples were assessed at an absorbance ratio of A 260/280 and A 260/230 , and the integrity of the RNA was confirmed by 1% agarose gel electrophoresis.

Library construction, Illumina sequencing and read assembly Fragment interruption, cDNA synthesis, addition of adapters, PCR amplification and RNA--Seq were performed by Beijing BioMarker Technologies (Beijing, China). The standard Illumina methods and protocols were adopted to prepare and sequence the cDNA libraries. NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB, USA) was used to enrich mRNA, and then the cDNA library was constructed using the NEBNext mRNA Library Prep Master Mix Set for Illumina (NEB, USA) and NEBNext Multiplex Oligos for Illumina (NEB, USA). The size of the library insert fragments was determined by 1.8% agarose gel electrophoresis, and the fragments were quantified using a Library Quantification Kit/Illumina GA Universal (Kapa, USA). Suitable fragments were selected as templates and sequenced on an Illumina HiSeq TM 2500 using paired-end technology. Three biological replicates were used to minimize sample differences.

In order to obtain clean and high-quality reads for sequence assembly, the raw reads were filtered by removing adaptor sequences, low-quality sequences (reads with ambiguous bases 'N') and reads with 10% > Q < 20% bases [27, 28] . The Trinity assemble program was used to assemble the clean reads into contigs, which covered more full-length transcripts through a broad range of expression levels [29] . The resultant contigs were added to transcripts based on paired-end information. The longest transcript from alternative splicing transcripts was selected as the unigene. These unigenes were combined to produce the final assembly and used for annotation.

To annotate unigenes, different sequences were searched by BLASTx against the NCBI nonredundant protein (nr) database and other databases, including Swiss-Prot protein database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Cluster of Orthologous Groups (COG) database. Gene Ontology (GO) annotations, including molecular functions, biological processes, and cellular components, were obtained using the Blast2GO program (https://www.blast2go.com/) [30, 31] . All searches were performed with an E-value < 10 −5 . Fragments per kilobase of transcript per million fragments mapped (FPKM) was calculated to represent the expression abundance of the unigenes [32] . FPKM may reflect the molar concentration of a transcript by normalizing for RNA length and for the total read number.

After normalizing genes expression levels, DEGs were obtained by pair-wise comparison of the four transcriptome libraries using IDEG6 software [33] . An FPKM fold change of 1.2 or 0.83 between two libraries was defined as the reference standard, with the Benjamini-Hochberg false discovery rate (FDR < 0.01) used to adjust the p-values.

In order to validate the results from our transcriptome sequencing analysis, the relative expression levels of 15 randomly selected genes were confirmed by reverse transcription quantitative PCR (RT-qPCR). Additionally, 9 genes with well reported previously were selected to further validate the genes of interest that might be involved in BmNPV resistance. All the Primers are listed in Table 1 . RT-qPCR reactions were prepared with the SYBR Premix Ex Taq TM Kit (TaKaRa), following the manufacturer's instruction. Reactions were carried out in Bio-Rad CFX96TM Real-Time System (Bio-Rad, USA). The thermal cycling profile consisted of an initial denaturation at 95°C for 30 s, 40 cycles at 95°C for 5 s, and 60°C for 30 s. All samples were performed in triplicate. Relative expression levels were calculated using the 2 -44Ct method following the protocol of Livak et al. [34] . In this study, B. mori ribosomal protein s3 (BmRPS3) gene was used as a reference gene. Statistical analysis was conducted using the SPSS software (IBM, www.ibm.com).

The LC 50 value was used to evaluate the resistant level of silkworm to BmNPV infection. The LC 50 value of A35 was approximately 26-fold greater than that of BC9 and over 500-fold greater than that of P50. The value of BC9 was 23-fold greater than that of P50 ( Table 2 ).

Transcriptome sequencing is an efficient technology for comparing gene expression levels in different samples, and in our study, it was used to search and analyze DEGs among P50, BC9 control and treatment groups. A total of four cDNA libraries were sequenced: P50-(P50 treated with sterile water), P50+ (P50 infected with BmNPV), BC9-(BC9 treated with sterile water), BC9+ (BC9 infected with BmNPV), with each group created in triplicate. After removing the adaptors and low quality sequences, 144,439,382 sequence reads were obtained ( Table 3 ). The GC content of each of the four libraries was approximately 50%, and CycleQ30% was greater than 89.91% for each library. Thus, the quality and accuracy of the sequencing data were sufficient for further analysis. Most of the reads matched silkworm genomic locations. All the unigenes matched previously described sequences with approximately 70% coverage. The length distribution of unigenes had similar patterns among the four libraries, suggesting that there was little bias in the construction of the four cDNA libraries (Fig 1) .

In order to annotate the unigenes, reference sequences were searched using BLASTX against the NCBI nr database (E-value < 10 −5 ). A total of 12,591 of 13,342 unigenes provided a BLAST result (S1 Table) . S2 Table shows the species with the closest match for each unigene. Most of Table 1 . Primers used in RT-qPCR for validation of DEGs.

Forward Primer Reverse primer

doi:10.1371/journal.pone.0155341.t001 

In order to determine the reliability of the transcriptome sequencing, the relative expression levels of 15 randomly selected genes were analyzed by RT-qPCR (Fig 2A) . The results were consistent with the transcriptome data. For example, the gene encoding peptidoglycan-recognition protein was down-regulated in both RNA-seq and RT-qPCR analyses, with a similar fold change. The lipase 1 gene was significantly up-regulated in the resistant strain BC9 after infection with BmNPV, which was also consistent with previously reported results [10] . Linear regression analysis of the correlation between RT-qPCR and RAN-seq ( Fig 2B) showed an R 2 (goodness of fit) value of 0.9169 and a corresponding slope of 1.5281, suggesting a strong positive correlation between RT-qPCR and transcriptome data. Therefore, the transcriptome data were satisfied for further analysis.

To further elucidate which DEGs had a potential role in antiviral response, a Venn diagram was constructed. A total of 285 DEGs were found to be differentially regulated when comparing P50+ and P50-, of which 122 were up-regulated and 163 were down-regulated. Similarly, 193 DEGs were found to be differentially regulated in the comparison of BC9+ vs. BC9-, with 56 genes up-regulated and 137 down-regulated. In addition, there were 154 DEGs differentially regulated in the comparison of BC9-and P50-, among which 78 genes were up-regulated and 76 genes were down-regulated (Fig 3, Table 4 , S3 Table) . There were 197, 119, 82 unique DEGs in P50+ vs. P50-, BC9+ vs. BC9-and BC9-vs. P50-, respectively. GO assignments were used to assign a functional classification to these DEGs. For cellular components, the number of unique DEGs fell into the macromolecular complex classification was distinct in BC9 and P50 following BmNPV infection. For molecular functions, the number of transporter activity related unique DEGs was distinct in BC9 and P50 following BmNPV infection. For biological progresses, the number of unique DEGs involved in metabolism processes and localization was distinct in BC9 and P50 following BmNPV infection (Fig 4) .

The comparisons of BC9+ vs. BC9-and P50+ vs. P50-identified many DEGs that might either be involved in silkworm defense against BmNPV or facilitate BmNPV infection. These genes could be divided into three categories: protein metabolism, cytoskeleton, and apoptosis.

Most of the DEGs (70%) associated with protein metabolism were down-regulated in BC9 following BmNPV infection. In contrast, 80% of the genes in P50 were up-regulated after BmNPV infection. Most of the DEGs (87.5%) associated with the cytoskeleton were up-regulated in P50 after BmNPV infection, but the number of up-regulated genes (50%) decreased in BC9. A total of 19 DEGs associated with apoptosis were identified following BmNPV infection. Ten of these DEGs were altered in BC9 after BmNPV infection. The other DEGs were all overexpressed in P50 following BmNPV infection (Table 5) .

Pathogen infection stimulated both cellular and humoral responses of insects [33] [34] [35] [36] . Genes participating in innate immunity pathways were identified and analyzed in regards to their potential role in BmNPV infection (Table 6 ). Thirty DEGs were identified, which could be classified into Toll pathway, Imd pathway, polyphenol oxidase (PPO) pathway, pattern recognition receptor, and antimicrobial peptide. As shown in Table 6 , 14 (47%) of these genes were downregulated and 9 (30%) were up-regulated in BC9 after infection with BmNPV, while 17 (57%) were down-regulated and 6 (20%) were up-regulated in P50 after BmNPV infection.

After removing the genetic background and strain specific immune stress response genes, a total of 22 DEGs were identified as possibly being involved in silkworm resistance to BmNPV (Fig 3, Table 7 ). For BC9-vs. P50-, 13 genes were up-regulated and the rest were down-regulated. However, the 22 genes were all down-regulated in BC9 following BmNPV infection ( Fig  5) . Some genes, including prostatic acid phosphatase, protease inhibitor 6, actin cytoskeleton- regulatory complex protein PAN1, and EF-hand domain-containing protein, were further analyzed by RT-qPCR (Fig 6) . After BmNPV infection, the expression levels of 4 genes were down-regulation in BC9 and A35 (resistant strain) (Fig 6) , which was consistent with the transcriptome data. Thus, we deduced that the 22 genes were possibly involved in resistance to BmNPV infection. Gene functions fell into the following categories: transport, virus replication, intracellular innate immunity, and apoptosis.

Despite the confirmation of an association between many genes and proteins and resistance to BmNPV, the molecular mechanism of antiviral activities was still unclear. Here, transcriptome sequencing was carried out to identify genes related to BmNPV-resistance in silkworm across the genome. By using the near-isogenic line BC9 (resistant strain) and the recurrent parent P50 (susceptible strain) to study silkworm antiviral mechanisms, some DEGs responding to BmNPV infection were successfully identified after comparing infected groups and controls in the two strains.

Based on the GO analysis, more DEGs were found to be involved in metabolic processes in BC9 (17.8%) than in P50 (9.2%) following BmNPV infection (Fig 4) , which was consistent with previous reports for Sogatella furcifera (Hemiptera: Delphacidae) and Campoletis sonorensis (Hymenoptera: Ichneumonidae) [35, 36] . Approximately twice as many DEGS related to transporter activity were identified in BC9 (12.2%) than in P50 (6.3%) (Fig 4) . Down-regulation of transporter related genes, such as lactase-phlorizin hydrolase (LPH), B(0, +)-type amino acid transporter 1 (BAT1), actin cytoskeleton-regulatory complex protein PAN1 (PAN1), MFS-type transporter (MFS), and otoferlin, could repress virus infection in host cells [37] [38] [39] [40] [41] [42] . Therefore, the increased number of these genes with altered expression levels in BC9 might be related to BmNPV infection. Moreover, the number of macromolecular complex genes in BC9 (14%) following BmNPV infection that showed an increase in expression compared with P50 (5.1%) were similar to previous reports [43, 44] .

Protein metabolism, cytoskeleton, and apoptosis may play an important role in host response to BmNPV infection Cellular and metabolic processes will be dramatically changed after viral infection [35] . In our study, several DEGs that participate in protein metabolism were found to be of interest. For example, solute carrier family 12 is involved in transporting extraordinarily diverse solutes [45] , and cystathionine gamma-lyase participates in hydrolysis of cystathionine [46] . Viruses may have to rely on cell proteins to accomplish replication in intercellular regions [47] , therefore, the downregulation of the genes involved in protein metabolism could inhibit the replication of BmNPV in the host cell. We speculated that the down-regulation of these genes affected virus replication. Cytoskeleton-dependent intracellular transport is an important strategy for transport of viral particles to different destinations [48] . In this study, some cytoskeleton related genes were found to be of interest, including actin cytoskeleton-regulatory complex protein PAN1 and actin-binding protein. These genes related to actin-coupled endocytosis could promote viral transport [49, 50] . We speculated that the down-regulation of the cytoskeleton genes might affect BmNPV transport.

Apoptosis plays a vital role in regulating cell response in Lepidopteran insects during viral infections, where larvae use selective apoptosis and subsequent sloughing of the infected cells in the midgut epithelium to resist virus infection [51, 52] . In this study, some genes related to activation of apoptosis were found to be of interest, including cytochrome c, inhibitor of caspase-activated DNase, amyloid precursor protein and B-cell lymphoma protein 2. The activity of caspase-activated DNase was blocked by Hepatitis C virus core at physiological levels, resulting in the inhibition of apoptotic cell death [53] . Amyloid precursor protein is a member of several signaling pathways that are involved in abnormal cell cycles, subsequently leading to apoptosis [54] . B-cell lymphoma protein 2 could bind to BH3 domains of various pro-apoptotic regulators to activate apoptosis [55] . The overexpression of cytochrome c in rabies virus showed a decreased pathogenicity in vitro and in vivo [56] . Based on their role in apoptosis activation, hosts need to increase the expression level of these genes to promote apoptosis when exposed to a virus; this supposition explains the up-regulation of genes involved in apoptosis in the transcriptome following BmNPV infection. Additionally, a significantly higher relative transcriptional level of cytochrome c was detected by RT-qPCR in BC9 and A35 (Fig  6) , which indicated the up-regulation of this gene could activate apoptosis to repress further BmNPV infection. However, some genes were down-regulated, including cysteine aspartic acid specific protease 9L, protein kinase A, apoptosis-inducing factor, apoptosis signal-regulating kinase 1, TNF-receptor-associated factor 6, TGF-beta-activated kinase 1, and p90 ribosomal S6 kinase. The inhibition of these genes may strongly impair viral infectivity and virus-induced apoptosis [57] [58] [59] [60] [61] [62] .

Although some differential expression of immune genes was observed, this might not be considered biologically important due to low expression levels; low expression levels were also found after BmCPV infection [21] . In this study, most immune related genes were down-regulated, with a few exceptions, such as the 2-fold up-regulation of toll-like protein and lysozyme. The expression of several other genes, including 18 wheeler, macrophage mannose receptor 1, and slit homolog 3 protein, were also up-regulated during the BmNPV infection. Unfortunately, the relationship between immune genes and BmNPV remains unclear and requires further study. We presumed that the low expression levels of immunity related genes may be associated with the disruption of the immune system by BmNPV, similar to the pathogenicity of human immunodeficiency virus (HIV).

Based on the Venn diagram, 22 genes of interest were identified, all of which were down-regulated in BC9 following BmNPV infection. These genes were grouped based on their functions, as reported in the literature. These groups, including transport, virus replication, intracellular innate immune, and apoptosis, may play an important role in the process of silkworm resistance to BmNPV. Several genes related to virus transport, including lactase-phlorizin hydrolase (LPH) [37] , B (0,+)-type amino acid transporter 1 (BAT1) [38, 46] , actin cytoskeleton-regulatory complex protein PAN1 (PAN1) [63] , and otoferlin [41, 42] were identified in this study. The expression level of these genes were all down-regulated in BC9 following BmNPV infection (Fig 5) . The resistant strain A35, a donor parent, was used to validate our results. The expression level of PAN1 was chosen for further testing by RT-qPCR. Expression levels of PAN1 were down-regulated in both BC9 and A35 following BmNPV infection (Fig 6) . The down-regulation of virus transport-related genes could inhibit the transmembrane and intracellular transport of BmNPV thereby preventing further infection. MFS-type transporter (MFS) [39, 40] , a transmembrane facilitator, was induced to increase intracellular monovalent ion concentrations, which led to lysis and cell death after HIV-1 infection. The down-regulation of this gene in silkworm could potentially block BmNPV infection.

Several genes related to virus replication were found to be of interest, including engrailed nuclear homeoprotein-regulated protein [64] , 1-acyl-sn-glycerol-3-phosphate acyltransferase alpha (ASGPA) [65, 66] , alanine aminotransferase 2 (ALT2) [67, 68] , U2AF domain protein (U2AF) [69, 70] , adenylate cyclase type 5 (ACT5) [71] , EF-hand domain-containing protein (EFHP) [72, 73] , prostatic acid phosphatase (PAP) [74] , and zinc ribbon domain protein (ZRDP) [75] . In this study, the expression level of these genes were all down-regulated in the transcriptome of BC9 following BmNPV infection (Fig 5) . In order to validate the results, RT-qPCR was conducted as described above. The expression levels of EFHP and PAP were lower following BmNPV infection in both BC9 and A35 (Fig 6) . The down-regulation of these genes could directly or indirectly participate in repressing BmNPV replication in host cells. Of the genes related to apoptosis, the expression level of pyruvate dehydrogenase kinase (PDK) was obviously inhibited after treatment with apoptosis-inducing agents [76] ; such down-regulation might activate apoptosis in response to BmNPV infection.

Protease inhibitor 6 is a member of the protein superfamily that contains TIL functional domain. Zhao et al. used genome sequences to demonstrate that the expression level of the TIL superfamily were down-regulated following BmNPV infection [77] , which was consistent with our results (Fig 6) . Furthermore, peptidoglycan-recognition protein 2, hemolin, facilitated trehalose transporter Tret1, integument esterase 2, defense protein precursor, and antimicrobial protein 6Tox precursor also showed differential expression following BmNPV infection. The relationship of these genes to BmNPV is an important area for further study.

Four DEGs were down-regulated in BC9 following BmNPV infection, including rasresponsive element-binding protein 1, GTPase-activating protein, trypsin alkaline A and peroxidase (Fig 5) . Previous studies revealed that these proteins play a role in virus infections. Mdv1-miR-M4 encoded by Marek's disease virus efficiently targeted the 3' untranslated regions of ras-responsive element-binding protein 1 (RREB1) [78] . TBC domain proteins belonging to a GTPase-activating protein were knocked out by double stranded RNA interference (RNAi), which led to a decrease in the level of transcripts of white spot syndrome virus genes [79] . Trypsin in the myocardium was able to trigger acute myocarditis following influenza A virus infection [80] . Overexpression and RNA silencing studies revealed that peroxidase was involved in the production of hepatitis C virus particles [81] . Moreover, RT-qPCR results indicated that the expression levels of all genes were down-regulated in BC9 and A35 following BmNPV infection (Fig 6) . Therefore, we speculated that the down-regulation of these genes might be involved resistance to BmNPV infection.

We hypothesized that the 22 DEGs discussed above played a role in the process of host response to BmNPV infection. The endocytosis process is triggered when the BmNPV nucleocapsid containing envelope binds to the cytomembrane. Vacuolar ATP synthase is activated by LPH to promote the fusion of the envelope and endosome thereby releasing the nucleocapsid into the cytoplasm. This process can be promoted by PAN1 and otoferlin. However, the transmembrane transport channel is an alternative pathway for virus to enter the cytoplasm, a process which can be facilitated by BAT1. The released nucleocapsid is transported into the nucleus with the help of the cytoskeleton (EFP). Once viral DNA is released into the nucleus, it will utilize host nucleotides to complete replication. In the final step of replication, viral DNA has to rely on host cell amino acids for assembly on the cytoskeleton (Table 5 ) [82] . In the cytoplasm, EFHP, ASGPA, ALT2, U2AF, ACT5 and ZRDP play an important role in facilitating virus replication, although the exact mechanism is still unclear. Moreover, transmembrane protein, MFS, is induced by BmNPV to increase cell volume, leading to lysis and cell death. We speculated that the down-regulation of these genes may affect the entry of virus into host cells and virus replication. The apoptosis process could also be triggered by PDK to inhibit BmNPV from further infecting other cells once BmNPV entered a host cell (Fig 7) .

Taken together, our results provide useful information on silkworm resistance to BmNPV infection.

Supporting Information S1 

The genetics and genomics of the silkworm, Bombyx mori

Silkworm genomics-progress and prospects

Genetic diversity among silkworm (Bombyx mori L., Lep., Bombycidae) germplasms revealed by microsatellites

Lipopolysaccharide elicits expression of immune-related genes in the silkworm, Bombyx mori

Gene expression profiling of resistant and susceptible Bombyx mori strains reveals nucleopolyhedrovirus-associated variations in host gene transcript levels

Screening of molecular markers for NPV resistance in Bombyx mori L. (Lep., Bombycidae)

Fluorescent differential display analysis of gene expression for NPV resistance in Bombyx mori L

Identification and characterization of an NPV infection-related gene Bmsop2 in Bombyx mori L

Antiviral activity of a serine protease from the digestive juice of Bombyx mori larvae against nucleopolyhedrovirus

A lipase isolated from the silkworm Bombyx mori shows antiviral activity against nucleopolyhedrovirus

In Vitro Antiviral Activity of an Alkaline Trypsin from the Digestive Juice of Bombyx Mori Larvae against Nucleopolyhedrovirus

Arginine kinase is highly expressed in a resistant strain of silkworm (Bombyx mori, Lepidoptera): Implication of its role in resistance to Bombyx mori nucleopolyhedrovirus

A Hypothetical Model of Crossing Bombyx mori Nucleopolyhedrovirus through Its Host Midgut Physical Barrier

Immune signaling pathways activated in response to different pathogenic micro-organisms in Bombyx mori

Expression of antimicrobial peptide genes in Bombyx mori gut modulated by oral bacterial infection and development

High-throughput transcriptome analysis of ISAV-infected Atlantic salmon Salmo salar unravels divergent immune responses associated to head-kidney, liver and gills tissues

Gut immunity in Lepidopteran insects

Hemocytes and hematopoiesis in the silkworm

Next-generation DNA sequencing techniques

De novo assembly and transcriptome characterization: novel insights into the natural resistance mechanisms of Microtus fortis against Schistosoma japonicum

Transcriptome analysis of Bombyx mori larval midgut during persistent and pathogenic cytoplasmic polyhedrosis virus infection

Genome-wide analysis of host gene expression in the silkworm cells infected with Bombyx mori nucleopolyhedrovirus

Transcriptome analysis of the brain of the silkworm Bombyx mori infected with Bombyx mori nucleopolyhedrovirus: A new insight into the molecular mechanism of enhanced locomotor activity induced by viral infection

Systemic and in vitro infection process of Bombyx mori nucleopolyhedrovirus

Detection of proliferation of Bombyx mori nucleopolyhedrovirus in its host by fluorescence quantitative PCR

Expression analysis of several antiviral related genes to BmNPV in different resistant strains of silkworm, Bombyx mori

De novo sequencing and a comprehensive analysis of purple sweet potato (Impomoea batatas L.) transcriptome

De novo assembly and characterization of root transcriptome using Illumina paired-end sequencing and development of cSSR markers in sweetpotato (Ipomoea batatas)

Full-length transcriptome assembly from RNA-Seq data without a reference genome

Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research

High-throughput functional annotation and data mining with the Blast2GO suite

Mapping and quantifying mammalian transcriptomes by RNA-Seq

IDEG6: a web tool for detection of differentially expressed genes in multiple tag sampling experiments

Analysis of relative gene expression data using real-time quantitative PCR and the 2(T)(-Delta Delta C) method

Transcriptome and Comparative Gene Expression Analysis of Sogatella furcifera (Horvath) in Response to Southern Rice Black-Streaked Dwarf Virus

Global Analysis of the Transcriptional Response of Whitefly to Tomato Yellow Leaf Curl China Virus Reveals the Relationship of Coevolved Adaptations

A Comparative Study on the Metabolism of Epimedium koreanum Nakai-Prenylated Flavonoids in Rats by an Intestinal Enzyme (Lactase Phlorizin Hydrolase) and Intestinal Flora

Differential influence of the 4F2 heavy chain and the protein related to b(0,+) amino acid transport on substrate affinity of the heteromeric b(0,+) amino acid transporter

A major superfamily of transmembrane facilitators that catalyse uniport, symport and antiport

Reduction of human immunodeficiency virus production and cytopathic effects by inhibitors of the Na+/K+/2Cl-cotransporter

Otoferlin is a calcium sensor that directly regulates SNARE-mediated membrane fusion

Ferlins: Regulators of Vesicle Fusion for Auditory Neurotransmission, Receptor Trafficking and Membrane Repair

Comparative proteome analysis of tracheal tissues in response to infectious bronchitis coronavirus, Newcastle disease virus, and avian influenza virus H9 subtype virus infection

Comparative study of the different activities of hepatitis B virus whole-X protein and HBx in hepatocarcinogenesis by proteomics and bioinformatics analysis

The ABCs of solute carriers: physiological, pathological and therapeutic implications of human membrane transport proteins-Introduction

Features of a Spatially Constrained Cystine Loop in the p10 FAST Protein Ectodomain Define a New Class of Viral Fusion Peptides

The cell cycle and virus infection

Intracellular Transport of Viruses and Their Components: Utilizing the Cytoskeleton and Membrane Highways

Functional analysis of Candida albicans genes whose Saccharomyces cerevisiae homologues are involved in endocytosis

Two novel members of the ABLIM protein family, ABLIM-2 and -3, associate with STARS and directly bind F-actin

Reaper homologue IBM1 in silkworm Bombyx mori induces apoptosis upon baculovirus infection

Insect defence: its impact on microbial control of insect pests

Antiapoptotic regulation by hepatitis C virus core protein through up-regulation of inhibitor of caspase-activated DNase

The cell cycle as a therapeutic target for Alzheimer's disease. Pharmacol Therapeut

Biophysical basis of the promiscuous binding of B-cell lymphoma protein 2 apoptotic repressor to BH3 ligands

Overexpression of cytochrome c by a recombinant rabies virus attenuates pathogenicity and enhances antiviral immunity

Influenza A Virus Induces Interleukin-27 through Cyclooxygenase-2 and Protein Kinase A Signaling

Porcine epidemic diarrhea virus induces caspase-independent apoptosis through activation of mitochondrial apoptosis-inducing factor

Mechanisms of minocycline-induced suppression of simian immunodeficiency virus encephalitis: inhibition of apoptosis signal-regulating kinase 1

Elevated suppressor of cytokine signaling-1 (SOCS-1): a mechanism for dysregulated osteoclastogenesis in HIV transgenic rats

Activation of NF-kappa B in virus-infected macrophages is dependent on mitochondrial oxidative stress and intracellular calcium: Downstream involvement of the kinases TGF-beta-activated kinase 1, mitogen-activated kinase/extracellular signalregulated kinase kinase 1, and I kappa B kinase

Phosphoproteomic Analysis of KSHV-Infected Cells Reveals Roles of ORF45-Activated RSK during Lytic Replication

Baculovirus gp64 envelope glycoprotein is sufficient to mediate pH-dependent membrane fusion

Requirement of the N-Terminal Activation Domain of Herpes Simplex Virus ICP4 for Viral Gene Expression

Enzymatic activities of the human AGPAT isoform 3 and isoform 5: localization of AGPAT5 to mitochondria

Phosphatidic Acid Produced by Phospholipase D Promotes RNA Replication of a Plant RNA Virus

cDNA cloning, genomic structure, chromosomal mapping, and functional expression of a novel human alanine aminotransferase

Activating transcription factor 4 mediates up-regulation of alanine aminotransferase 2 gene expression under metabolic stress

Genomic functions of U2AF in constitutive and regulated splicing

A suboptimal src 3' splice site is necessary for efficient replication of Rous sarcoma virus

Herpesvirus quiescence in neuronal cells IV: virus activation induced by pituitary adenylate cyclase-activating polypeptide (PACAP) involves the protein kinase A pathway

EFHC1 interacts with microtubules to regulate cell division and cortical development

Selective reovirus infection of murine hepatocarcinoma cells during cell division. A model of viral liver infection

Naturally occurring fragments from two distinct regions of the prostatic acid phosphatase form amyloidogenic enhancers of HIV infection

Zinc Ribbon Domain-Containing 1) Is a Host Cellular Factor That Influences HIV-1 Replication and Disease Progression

Pyruvate dehydrogenase kinase 4 (PDK4) could be involved in a regulatory role in apoptosis and a link between apoptosis and insulin resistance

Genome-Wide Identification and Immune Response Analysis of Serine Protease Inhibitor Genes in the Silkworm, Bombyx mori

Marek's disease virus microRNA designated Mdv1-pre-miR-M4 targets both cellular and viral genes

PmTBC1D20, a Rab GTPase-activating protein from the black tiger shrimp, Penaeus monodon, is involved in white spot syndrome virus infection

Ectopic trypsin in the myocardium promotes dilated cardiomyopathy after influenza A virus infection

Quantitative Proteomics Identifies the Membrane-Associated Peroxidase GPx8 as a Cellular Substrate of the Hepatitis C Virus NS3-4A Protease

Baculovirus infectivity and the actin cytoskeleton