Transcriptome Analysis of Bombyx mori Larval Midgut during Persistent and Pathogenic Cytoplasmic Polyhedrosis Virus Infection

Many insects can be persistently infected with viruses but do not show any obvious adverse effects with respect to physiology, development or reproduction. Here, Bombyx mori strain Daizo, persistently infected with cytoplasmic polyhedrosis virus (BmCPV), was used to study the host’s transcriptional response after pathogenic infection with the same virus in midgut tissue of larvae persistently and pathogenically infected as 2nd and 4th instars. Next generation sequencing revealed that from 13,769 expressed genes, 167 were upregulated and 141 downregulated in both larval instars following pathogenic infection. Several genes that could possibly be involved in B. mori immune response against BmCPV or that may be induced by the virus in order to increase infectivity were identified, whereas classification of differentially expressed transcripts (confirmed by qRT-PCR) resulted in gene categories related to physical barriers, immune responses, proteolytic / metabolic enzymes, heat-shock proteins, hormonal signaling and uncharacterized proteins. Comparison of our data with the available literature (pathogenic infection of persistently vs. non-persistently infected larvae) unveiled various similarities of response in both cases, which suggests that pre-existing persistent infection does not affect in a major way the transcriptome response against pathogenic infection. To investigate the possible host’s RNAi response against BmCPV challenge, the differential expression of RNAi-related genes and the accumulation of viral small RNAs (vsRNAs) were studied. During pathogenic infection, siRNA-like traces like the 2-fold up-regulation of the core RNAi genes Ago-2 and Dcr-2 as well as a peak of 20 nt small RNAs were observed. Interestingly, vsRNAs of the same size were detected at lower rates in persistently infected larvae. Collectively, our data provide an initial assessment of the relative significance of persistent infection of silkworm larvae on the host response following pathogenic infection with CPV, while they also highlight the relative importance of RNAi as an antiviral mechanism.


Introduction
Interest in viruses originates mainly from their role as pathogens that can cause significant disease in hosts. Regarding pathogenic infections in insects, virus research focuses on the protection of beneficial insects such as the honeybee, Apis mellifera, and the silkworm, Bombyx mori, against virus infection [1,2] as well as on the improvement of the use of viruses to control agricultural pests [3]. In addition, many research activities focus on the transmission by mosquitoes of arbovirus infections that cause disease in humans and livestock [4,5].
With the advent of next generation sequencing technologies, however, it has been realized that many insect species can be persistently infected with viruses that do not cause an obvious adverse symptom or abnormal phenotype [6,7]. In contrast to acute infections, which are limited in time, resulting either in clearance of the virus or death of the host, persistent infections can be maintained in the hosts for very long periods and transmitted to the offspring [8]. With time, mutualistic relationships can even be developed between host and virus and, in some cases, it has been documented that viral genes and mechanisms are adopted by the host for its own benefit, for instance in the immune response [9]. Of great interest are also arbovirus infections in mosquito vectors, which can be considered as persistent infections where both mosquito vector and virus need to survive for a sufficient time to allow transmission to the vertebrate link in the infection cycle [10]. Also honeybee populations can be persistently infected with multiple viruses without affecting the health of the colony [11]. However, bee viruses are often mentioned as one of the multifactorial causes of the decline of honeybee colonies [12] and it can therefore be hypothesized that persistent infections could turn pathogenic in conditions of increased stress. Because of their prevalence, persistent infections deserve more attention with respect to alterations in the physiology of the host, such as the immune response, and to conditions that could turn apparently harmless persistent infections into disease.
B. mori cytoplasmic polyhedrosis virus (BmCPV), a reovirus characterized by a segmented dsRNA genome, is a major pathogen of the silkworm [13]. In contrast to baculovirus infections, which can spread efficiently throughout the body and are very virulent to the host, CPV infections are mostly or exclusively limited to the midgut tissue, cause less damage and therefore have a propensity to become persistent [14]. During our research with BmCPV, it was noticed that our laboratory strain of silkworm, Daizo, was persistently infected with BmCPV. Despite the infection, however, animals appeared healthy, showed normal growth and metamorphosis, and produced batches of fertilized eggs of typical size. Because of this observation, we considered the persistently infected silkworm strain as an interesting experimental model to investigate the (alteration of the) transcriptional response after feeding of a high dose of BmCPV polyhedra that caused clear pathogenic effects. Several studies concerning the transcriptional responses to BmCPV infection can be found in literature; however, these responses were recorded in larvae for which no persistent infection was reported [15][16][17][18]. Another study was focused on the miRNAs' differential expression during BmCPV infection [19].
In our study, deep sequencing (Illumina) technology was applied to obtain an initial assessment of the transcriptional response in silkworm larvae that were persistently or pathogenically infected with BmCPV. Instead of having different biological replicates at a single developmental stage, it was decided in this explorative study to compare the transcriptional responses between persistently and pathogenically infected animals at two different developmental stages (2 nd and 4 th larval instars) and to focus on genes that become differentially expressed irrespective of the developmental stage. The identification of differentially expressed genes was subsequently validated by qRT-PCR experiments on samples obtained from 2 nd instar larvae. Despite the limitation of the deep sequencing analysis (data from four unique samples, corresponding to persistent or pathogenic infection in the 2 nd or 4 th instar), it is believed that valuable preliminary data were obtained that will stimulate further research. It is noted that a similar format of experimental design was applied to the study of differentially expressed microRNAs (miRNAs) after BmCPV infection [19].
Our study establishes a considerable overlap in the transcriptional response to BmCPV during a pathogenic infection of persistently infected silkworm larvae with the transcriptional response documented in the above mentioned studies where no persistent infection was reported. In addition, we focused on the involvement of the RNA interference (RNAi) machinery during BmCPV infection, which has not received much attention in previous BmCPV-related studies. Although RNAi has been considered as the most important antiviral response in Drosophila and mosquitoes [20][21][22], in insects of other groups, such as the honeybee or the silkworm, its involvement in antiviral defense remains largely unknown. Attention is therefore paid to the transcriptional response of RNAi machinery genes during a pathogenic BmCPV infection as well as to the detection of virus-derived small RNAs (vsRNAs).

Silkworm rearing and infection with BmCPV
The larvae of B. mori, Daizo strain, were reared on artificial diet (Yakuruto, Tokyo, Japan) at 25°C under a photoperiod of 12 h light and 12 h dark. For pathogenic infection with BmCPV, aliquots of artificial diet were treated with 50 μl of a concentrated solution of polyhedra of BmCPV (3.25 x 10 7 polyhedra/ml) and fed to 2 nd instar (N = 30) or 4 th instar larvae (N = 30), 1-2 days after molt. Successful pathogenic infection was indicated by retardation of growth and confirmed by detection of polyhedra by means of optical microscopy as well as by genespecific RT-PCR for polyhedrin (see Results and Discussion section). Midgut and body wall tissue was collected from pathogenically infected larvae at 10-20 days after feeding of BmCPV polyhedra. The same tissues were also collected from control (persistently infected) larvae at a similar stage.

RNA extraction and reverse-transcription PCR (RT-PCR)
Tissues were homogenized in TRI Reagent (Sigma, Saint Louis, MO) and total RNA was extracted according to the manufacturer's protocol. The quantity of extracted RNA was assessed with a NanoDrop 1000 Spectrophotometer (Thermo Scientific, Waltham, MA) and/or by electrophoresis on 1% (w/v) agarose gels. RNA for specific detection of BmCPV polyhedrin was first mixed with DMSO at 1:1 ratio and heated at 50°C for 45 minutes for denaturation of dsRNA BmCPV genome. One microgram of total RNA was used as template for first-strand complementary DNA (cDNA) synthesis as performed by a Revert Aid reverse transcriptase (Thermo Scientific). Random hexamers (Thermo Scientific) were utilized as primers for cDNA templates of quantitative RT-PCR, whereas BmCPV segment 10-specific primer 5'-AGGAT-CATGGCAGACGTAGC-3' was used to synthesize cDNA templates for polyhedrin-specific RT-PCR. Viral polyhedrin was detected using the same RT primer as forward primer and 5'-TAGGCGTTCGGCGAAATGT-3' as reverse primer, under the following cycling conditions: 94°C for 30 s, 50°C for 30 s and 72°C for 45 s (40 cycles).

Deep sequencing of RNA samples from midguts of pathogenically and persistently infected larvae
The four RNA samples that were prepared for deep sequencing analysis were derived from midgut tissue of persistently ("control") and pathogenically infected 2 nd instar larvae (2c and 2inf), as well as of persistently ("control") and pathogenically infected 4 th instar larvae (4c and 4inf). RNA quantifications were performed using Qubit fluorometry (Life technologies, Carlsbad, CA).
For each sample, an Illumina mRNA sequencing library was made from 100 ng of total RNA using the TruSeq Stranded mRNA Sample Prep Kit (Illumina, San Diego, CA), whereas approximately 500 ng of total RNA was used to create a small RNA library using the TruSeq Small RNA Sample Preparation Kit (Illumina). The 4 mRNA and the 4 small RNA libraries were each equimolarly pooled and sequenced in one lane of an Illumina HiSeq 2000 flowcell, generating 1 x 50 bp reads. After sequencing, the data was demultiplexed using the sample specific nucleotide barcodes. On average, 30 x 10 6 mRNAs were generated. The mRNA and small RNA differential expression analysis was performed using CLC bio (Qiagen, Venlo, The Netherlands). All reads were trimmed for Illumina adapter sequences.
The mRNA reads were mapped to the Bombyx annotated genome (Kaikobase: http://sgp. dna.affrc.go.jp/pubdata/genomicsequences.html) [23]. For the 2inf, 2c, 4inf and 4c mRNA samples, the percentage of mapped reads was 24%, 37%, 29% and 36%, respectively. The lower percentages of reads that were mapped in the infected samples likely reflect the predominance of transcripts that could be mapped to the viral genome, which illustrates the severity of the infection in those samples. The number of mRNA reads that mapped to a transcript was divided by the transcript length and normalized per sample by the number of mapped reads to calculate the RPKM (Reads Per Kilobase per Million mapped reads) expression values for each gene, thus making possible the direct comparison of differential expression among samples [24,25]. Because of the lack of biological repeats, we did not use T statistics to identify significant differentially expressed genes. Instead, we filtered the RNA-seq data according to the following criteria: (i) total gene reads count resulting from all 4 libraries should be more than 10, and (ii) selected transcripts should present at least a 2-fold RPKM change (|log2Ratio| 1) between pathogenically and persistently infected in both library pairs (2inf/2c, 4inf/4c) towards the same direction, i.e. up-or down-regulation [26,27]. These criteria therefore ensure, first, that all relevant genes including those with low, but significant, expression level are evaluated, and, second, that a significant change in expression occurs that is independent of the developmental stage of the infected larvae.
Regarding small RNAs, after trimming, more than 95% of the original number of reads was used for small RNA tag counting. Approximately 15% of the tags had a count lower than 5 and these were discarded. The remaining tags were counted and annotated using the Bombyx miRNA miRBASE 19 database. The small RNAs that did not map against the Bombyx miRNA miRBASE 19 database were subsequently mapped against the BmCPV genome to find vsRNAs originating from the virus. Small RNA length distribution graphs were produced in R [28] by use of the viRome package (http://www.ark-genomics.org/bioinformatics/virome). The sequences from the mRNA and small RNA libraries of all 4 samples were submitted to the European Nucleotide Archive (Accession Number PRJEB7502).

Gene Ontology (GO) analysis
A list of sequences that were identified as differentially expressed between persistently and pathogenically infected larvae (see Results and Discussion section) were analyzed for GO annotation [http://www.geneontology.org [29]]. Sequences were identified based on the annotated gene set at Kaikobase and genes resulting from the homology searches were used as input in DAVID [http://david.abcc.ncifcrf.gov [30]].

Real-time quantitative RT-PCR
Real-time qRT-PCR was performed in a Mx3000P QPCR System (Agilent Technologies, Santa Clara, CA) equipped with MxPro QPCR Software (Agilent Technologies) using KAPA SYBR FAST qPCR Kit (Kapa Biosystems, Wilmington, MA), gene-specific primers at a final concentration of 0.55 μM each and 12.5 ng of midgut cDNA template. Relative transcript levels were normalized to the expression level of the cellular actin gene. PCR cycling started with initial activation of KAPA SYBR FAST Master Mix polymerase at 95°C for 3 min, followed by 40 cycles of 95°C for 5 s, 59 or 60°C for 30 s and 72°C for 5 s. Forward and reverse primers to detect specific transcripts (S1 Table) were designed using Primer Express 2.0 software (Applied Biosystems/ Life Technologies). Relative expression levels of target gene (X) were calculated in relation to the transcription levels of the actin reference gene (R), as 2 -ΔCt , where ΔC t = C t X -C t R .

Persistent (non-pathogenic) infection of silkworm with BmCPV
To verify infection of silkworm larvae with BmCPV, we performed RT-PCR with gene-specific primers for the viral polyhedrin gene. The samples tested corresponded to RNA collected from body wall, midgut, midgut content and feces. For pathogenic infections (2inf and 4inf samples), polyhedrin product was detected at high levels in all samples analyzed (midgut tissue, midgut content, feces and body wall; Fig. 1). However, rather high amounts of polyhedrin RNA could also be detected in the midgut content of the control (not pathogenically infected) 2 nd instar larvae (2c sample) and low levels were also present in the midgut tissue (both 2 nd and 4 th instar stages; Fig. 1). In addition, cDNA samples derived from the eggs of the same strain as the larvae (Daizo) were also positive for polyhedrin RNA, in contrast to the eggs of the P50 strain, which gave no PCR product (Fig. 1). PCR products obtained from control midguts (4 th instar) and Daizo eggs were validated by sequencing to correspond to BmCPV polyhedrin. The detection of polyhedrin RNA in midgut tissue and eggs of animals of the Daizo strain in our colony is indicative of a persistent state of infection by BmCPV.
Thus, larvae of the Daizo strain that were pathogenically infected in our study with BmCPV were already persistently infected with the same virus at low levels. This persistent state of infection was not visibly correlated with pathogenic effects, as larvae apparently grew, molted and metamorphosed normally. It should also be mentioned that polyhedra were never detected in any control animal under normal rearing following microscope observation. Furthermore, RT-PCR experiments revealed the presence of polyhedrin RNA mainly in the midgut content of control larvae and only traces in the midgut tissue itself (Fig. 1), which suggests that in persistently infected larvae, the infection of the midgut tissue (epithelium and associated muscle and connective tissue) is efficiently cleared, thus preventing the occurrence of pathogenic defects. This contrasts to the situation in pathogenically infected larvae, where high amounts of polyhedrin RNA were detected in the midgut tissue ( Fig. 1), being indicative of high levels of viral replication/production (see further below). Finally, virus-derived small RNAs (vsRNAs) were detected in both persistently and pathogenically infected larvae ( Fig. 2; see also Detection and preliminary analysis of vsRNAs paragraph); however, the detection levels in the control larvae were minimal, further corroborating the low levels of viral RNA in persistently infected tissues.

Oral infection of silkworm larvae with BmCPV polyhedra
Persistently infected larvae of Daizo strain were pathogenically infected through the oral route with 5.5 x 10 4 polyhedra per larva. At the time of infection, the larvae were at the beginning of the 2 nd or the 4 th instar developmental stage. The infection was verified macroscopically by observation of the infected larvae (S1 Fig.), as well as microscopically by detection of viral polyhedra in the larval midgut tissue (S2 Fig.). As shown for a representative group of persistently and pathogenically infected larvae in S1 Fig., the sizes of the pathogenically infected larvae varied notably among several individuals. This phenomenon was observed in both 2 nd and 4 th instar infected larvae. However, polyhedra could be detected in larvae of all sizes, indicating that the impact of the infection on body size did not correlate with the presence of infection.
According to the literature, BmCPV targets epithelial midgut cells and is not expected to expand all over the larval body. Nevertheless, following observation of several larval organs under the microscope, we have detected polyhedra also in the interior of the body wall, as well as in the hemolymph (S2 Fig.).

Transcriptome analysis of larval midgut samples
RNA-seq was used to analyze midgut samples from silkworm larvae infected in the 2 nd and 4 th instar stage, in order to detect differentially expressed transcripts during persistent and pathogenic BmCPV infection. Deep sequencing was carried out on samples from two different developmental stages to identify transcripts that are differentially expressed irrespective of the developmental stage, instead of having biological replicates at the same stage of development. Transcriptome analysis resulted in the identification of 13,769 expressed genes over all 4 RNA samples (S1 Dataset). Filtering of RNA-seq data according to the criteria outlined in Materials and Methods section resulted in 167 up-regulated and 141 down-regulated transcripts in both instars ( Fig. 3; S2 and S3 Datasets).
According to GO analysis of the 308 differentially expressed transcripts, 201 were found to have a GO annotation and were initially categorized in the three main GO functional groups, i.e. biological process, molecular function and cellular component. Genes belonging to each group were further classified at level 2 for each functional group (Fig. 4).
Gene transcripts found to have at least one"biological process" annotation (426), were further classified in 17 subgroups, with the most important of them belonging to metabolic process (21.1%), cellular process (18.3%), single-organism process (15%), response to stimulus (10.6%) and multicellular organismal process (10.1%) ( Fig. 4; S3 Fig.). Of note is that the category "immune system process" was poorly represented. Gene transcripts with annotations of molecular function (265) were further classified in 10 categories, among which catalytic activity Differentially expressed genes and their possible connection to a BmCPV-specific antiviral response Several genes possibly implicated in silkworm's response against BmCPV infection were found by deep sequencing to be highly differentially expressed between persistently and pathogenically infected samples (S2 and S3 Datasets). Some of the most interesting genes were further validated by qRT-PCR ( Fig. 5; S4 Fig.). These genes fall in several categories, i.e. physical barriers, immune responses, proteolytic enzymes, heat-shock proteins (Table 1), metabolic enzymes, hormonal signaling and uncharacterized, as outlined in detail in Table 2. Together, they constitute a complex response to BmCPV infection in the silkworm larvae.
In the literature, several reports exist with respect to the transcriptome response to BmCPV infection using larvae without persistent virus infections [15][16][17][18]. Table 3 presents a list of all differentially expressed genes obtained from previously published studies and from our study. This list contains only genes of which the differential expression was confirmed by qRT-PCR. In accordance with all four studies, down-regulation of expression of several digestive enzymes was noted; an induction of a small set of immune response genes, of which insulin-binding protein 2 (ibp2) gene was identified as significantly induced in three independent studies ( [16,18]; our study); an induction of the apoptosis pathway or repression of apoptosis inhibitors in four studies ( [15][16][17]; our study) was shown, as well as a role for remodeling of the midgut epithelium during infection (Table 3). For the latter process, it has been speculated that the induction of the cuticle protein-like gene CPH43, as observed in our study, plays a role. In two studies,  Indicated is the fold change in expression of selected genes between persistent and pathogenic infection as obtained by qRT-PCR on 2 nd instar samples (qPCR-2) or by deep sequencing on both 2 nd and 4 th instar samples (DS-2 and DS-4, respectively). Selected genes belong to the group with highest difference in expression between persistent and pathogenic infection. Please note that fold changes are expressed as log2 values (a 2-fold up-or downregulation corresponds to a log2 value of 1 or -1, respectively). See Table 2  Immune responses

IBP2
highly homologous to IGFBP7 in vertebrates (tumor-suppressive function leading to apoptosis); involved in immune and endocrine responses [32][33][34][35][36][37] insulin signaling pathway is reported to be involved in the defense against pathogens in C. elegans [38] vertebrate insulin activates the extracellular signal-regulated kinase (ERK) in the mosquito gut; has an antiviral role in Drosophila cells and insect gut epithelium [39] PGRP LB-like behaves as an amidase; hydrolyzes Grambacteria peptidoglycan and is activated by them; indirectly acts as negative regulator of the IMD pathway, thus balancing homeostasis and immune response activation [40,41] Proteolytic enzymes Caspase-8-like death receptor-associated initiator caspase-8 is expected to be activated upon reovirus infection in the case of reovirus-induced apoptosis [42][43][44][45] CPB up-regulated in the midgut of Anopheles gambiae mosquitoes infected with Plasmodium falciparum parasite, significantly reduced growth and development of rodent parasite P. berghei in mosquitos fed on infected mice immunized against CPB [46] Trypsin-like protease precursor alkaline protease; active in the silkworm midgut alkaline environment; participates in the hydrolysis of incoming food and probably also of viral polyhedra, thus mediating release of occluded viruses and infection of midgut columnar epithelial cells [47] trypsin-like protease is down-regulated in BmDNV-Z and BmNPV-infected silkworm strains susceptible to the respective viruses, up-regulated in BmNPV-infected silkworm strain [47,48] facilitates DENV-2 virus infection in Aedes aegypti [49] Zinc carboxypeptidase A1-like (CPA1) Metallocarboxypeptidase ACI, a zinc carboxypeptidase A1 (CPA1) inhibitor, is expressed by Ascaris (human intestinal parasite) to enhance its survival during infection [50] Heat-shock proteins HSPs; HSCs molecular chaperones in various cellular processes; danger signals thus activating host immune response [51,52] HSPs and HSCs are necessary for efficient BmNPV proliferation in Bombyx cells as well as for PCV2 virus expansion in porcine cells [53][54][55] sHSPs are induced in BmCPV-infected larval midguts [18] Hsp25 has antiviral role in reovirus-infected murine cells [56]; HSC70 is up-regulated in BmCPV-infected silkworms (72hpi) [18]; Hsc70t and HSP105 are induced by reovirus infection in murine cells [57] Metabolic enzymes DCXR converts L-xylulose in xylitol (carbohydrate metabolism); reduces the highly reactive α-dicarbonyl compounds (DCs) with endogenous/ xenobiotic origin (detoxifying enzyme) [58] DHS-21 (DCXR ortholog) is essential for normal life-span and reproduction of C. elegans [59] Esterase FE4-like esterases hydrolyze and inactivate insecticides (insecticide resistance), but also metabolize pathogen-secreted toxic compounds (host response) [60]; as counter-mechanism, pathogens may synthesize inhibitors against such esterase function to cause down-regulation of esterase genes Hormonal signaling NPLP4E neuropeptides: signaling molecules playing key roles in insects due to their involvement in developmental, reproductive, metabolic and behavioral processes; expressed in several tissues (brain, epidermis, ovary, prothoracic gland); possibly involved in molting regulation [61] NPLP4E maybe carries out additional functions in the midgut (here was first detected to be expressed) [31,62] Uncharacterized LOC101735726 possible role in the physiological process of DNA replication and maintenance (S5 Fig.) which may be hindered due to BmCPV presence in the midgut The list represents genes with high differential expression levels during pathogenic infection detected by deep sequencing analysis that were also validated by qRT-PCR.  Genes SH (p50) [17] MA (p50) [18] DS (4008) [15] DS (4008) [16] DS (p50) [ [15,16] also emphasized the possible role of the calcium-signaling pathway to induce apoptosis in infected cells, but the calcium signaling pathway genes were not among the most differentially genes that were analyzed by qRT-PCR in the other studies. The differences in transcriptional response among the four studies can be explained by the use of different silkworm strains, stage of infection and time of collection of infected samples. Gao et al. [16] compared the transcriptional response to BmCPV between a relatively resistant Genes SH (p50) [17] MA (p50) [18] DS (4008) [15] DS (4008) [16] DS (p50) [16] D (DaiSzo) our study (p50) and a susceptible (4008) silkworm strain and found important differences in the response to virus infection. Of interest, it was observed that ibp2 gene was induced at high levels only in the resistant strain. Therefore it was considered as a gene involved in resistance against BmCPV. The observation of induction of ibp2 during pathogenic infection in our study suggests the existence of an active antiviral resistance mechanism during persistent infection with BmCPV.
Transcriptional response of RNAi-related genes to infection with BmCPV in B. mori As outlined in the introduction, a second goal in our study was to assess the RNAi response during a persistent and pathogenic infection with BmCPV since this response did not receive much attention in previous studies. Therefore, the transcriptional response of RNAi-related genes in the host tissues was examined, while the production of vsRNAs was also investigated. It is well known that in Drosophila and mosquitoes, the siRNA pathway (and piRNA pathway in mosquitoes) can play an active role in the defense against exogenous invading singlestranded and double-stranded RNAs [63]. Due to the fact that the BmCPV genome consists of 10 segments of dsRNA, it was speculated that it could act as a target for the RNAi mechanism during the establishment of pathogenic BmCPV infection in Bombyx larvae. Therefore, it would be interesting to investigate whether pathogenic infection could up-regulate the expression of RNAi-related genes as an antiviral response mechanism during pathogenic infection.
Following application of the BLASTP algorithm in the NCBI database on the B. mori genome, as well as use of the BLAST tool in the www.silkdb.org database, homologs (annotated or not) of several RNAi-related genes originally implicated in the RNAi process in other organisms were identified in the silkworm's genome (for a discussion of different categories of RNAi-related genes, see [64]). Deep sequencing data were then analyzed to pinpoint alterations of the expression of RNAi-related genes following a pathogenic infection with BmCPV. The expression of several core RNAi genes and additional RNAi-related factors (Table 4) was further confirmed by qRT-PCR in 2 nd instar midgut cDNA samples ( Fig. 6; S6 Fig.).
Three genes of the Argonaute family (i.e. Ago1, Ago2 and Ago3 that represent the miRNA, siRNA and piRNA pathway, respectively), as well as Dcr2 from the siRNA pathway, are considered as core RNAi genes in B. mori, since inhibition of their expression by use of the "RNAi-ofthe RNAi" method caused a significant impediment of the silencing potency in the silkwormderived Bm5 cells [65]. In addition, several other core RNAi genes, such as Aubergine (SIWI; piRNA pathway) and Loquacious (dsRNA binding protein; miRNA pathway), were also analyzed in the deep sequencing database and by qRT-PCR (Table 4; Fig. 6; S6 Fig.).
Regarding the siRNA pathway, the genes Ago2 and Dcr2, which are well known to be responsible for the defense against exogenous dsRNAs in Drosophila, showed a~2-fold up-regulation in expression during the pathogenic infection (Table 4; Fig. 6; S6 Fig.). The piRNA pathway gene Aubergine was found to be 1.5-fold up-regulated in the 4 th instar animals after feeding with polyhedra; however, this result was not confirmed for the 2 nd instar stage. Also for Loquacious (miRNA pathway), no important differential expression between the two types of infection was observed (Table 4; Fig. 6; S6 Fig.).
Apart from the core RNAi components, additional RNAi-related genes were selected [64]. Tudor staphylococcus / micrococcal nuclease (Tudor-SN) and eukaryotic initiation factor-4E1 (eIF-4E1) were recently shown to constitute core factors for the formation of stress granules (SGs) in silkworm BmN4 cells, while they possibly also interact with Ago1 and Ago2 proteins of the miRNA and the siRNA pathway, respectively [66]. Tudor-SN expression was not notably altered after a pathogenic infection with BmCPV (Table 5; Fig. 6; S6 Fig.). Interestingly, however, the expression level of eIF-4E1 was up-regulated during the pathogenic infection (from 2-to 8-fold, depending on the technique used and the sample; Table 5; Fig. 6; S6 Fig.). The increase of the relative mRNA level could be due to the general stress state caused by the infection, or due to a special role of eIF-4E1 in the midgut response.   Table 5 and [64] for further explanation on gene identity and function.
doi:10.1371/journal.pone.0121447.g006 Heat shock protein 90 (HSP90) is an important factor that lately has emerged as a key piRNA pathway regulator. HSP90 was shown to take part in the piRNA pathway in BmN4 cells, more specifically in the recruitment of precursor piRNA molecules by the PIWI proteins [67,68]. Although deep sequencing analysis showed a >2-fold up-regulation of HSP90 in the 4 th instar library of pathogenically infected midgut tissue, this was not confirmed in the respective 2 nd instar sample (Table 5; Fig. 6; S6 Fig.). These data indicate that the role of HSP90 is relatively minor in this type of infection, in sharp contrast to other HSPs as shown in Tables 1-3.
Regarding genes involved in dsRNA uptake, CG4572 and HPS4 only showed clear differential expression between pathogenic and persistent infection in the 4 th instar animals. Also for the moderately expressed scavenger receptor-C (SR-C), no clear differential expression was observed (Table 5). Ratios of RPKM values from pathogenically versus persistently infected midgut tissue of 2 nd and 4 th instar larvae obtained by deep sequencing analysis (2c, 2inf, 4c and 4inf samples) are presented for selected genes with a role in the RNAi mechanism. Listed are intracellular auxiliary factors, dsRNA uptake genes, antiviral RNAi genes, nucleases and unclassified factors [64]. Genes presenting higher than 1.5-fold up-or down-regulation are marked with bold letters. Abbreviation: nd: not detected. doi:10.1371/journal.pone.0121447.t005

Alteration in immune gene expression after infection with BmCPV
Genes belonging to innate immunity pathways were identified and analyzed regarding their possible differential expression upon infection (Table 6). Although differences were observed, these could often not be considered as biologically important because of low expression levels (RPKMs). However, a few exceptions were noted, such as the 2-fold down-regulation of Toll6 and of one Attacin1 homologue, as well as the 2-fold up-regulation of beta-1,3-glucan recognition protein 2 gene in pathogenically infected larvae. The expression of cecropins (A, B and E) also tended to be up-regulated during the pathogenic infection. On the other hand, Toll9-1, which was previously found to be down-regulated by exogenous dsRNA application [69], was not differentially expressed in any of the two pathogenically infected midgut samples.

Detection and preliminary analysis of viral small RNAs (vsRNAs)
Deep sequencing was also used for the analysis of the small RNAs in samples of persistently and pathogenically infected larvae. Filtering of the small RNAs for reads that map to the BmCPV genome resulted in a total of 4,487,417 reads for all 4 samples. When the number of reads of the vsRNAs was plotted against their length, a clear peak of 20 nt vsRNAs equally distributed between both genomic dsRNA strands appeared in all 4 samples. As expected, vsRNA reads were highly abundant in pathogenically infected midguts, while for persistently infected samples their abundance was notably low. The peak of 20 nt corresponds to 749,372 and 2,079,497 reads for the 2inf and 4inf samples, respectively, while for the control samples 2c and 4c the small RNA counts were 324 and 481, respectively (Fig. 2). These data clearly indicate that the RNAi machinery is responding to BmCPV infection and that the abundance of produced vsRNAs is correlated with the severity of the infection. Further analysis of the vsRNA data such as the mapping of the vsRNAs to viral dsRNA genome segments to detect hot-spots of small RNA production and to identify differences in vsRNAs between persistently and pathogenically infected animals is currently being carried out (manuscript in preparation). The observation that the vsRNAs mapped equally to sense and antisense strands of the dsRNA genome (Fig. 2) strongly indicates that the dsRNA genome segments, rather than structured dsRNA regions of viral (sense) mRNAs, are the source for production of vsRNAs by Dicer enzymes. Further studies are required to establish the functionality of the vsRNAs and to investigate whether differences in activity exist between vsRNAs produced in persistently versus pathogenically infected animals as well as between different regions of the viral dsRNA genome (cold-spots versus hot-spots).

Conclusions
In this work, the discovery of persistent BmCPV infection of our silkworm laboratory colony provided an opportunity to compare the transcriptional response to pathogenic infection with that occurring in non-persistently infected larvae, as described in the literature. Our conclusions can be summarized as follows: 1. The transcriptional response against pathogenic BmCPV infection is complex and suggests the involvement of several mechanisms, including RNAi.
2. Pre-existing persistent infection does not profoundly affect the antiviral response against pathogenic infection with the same virus, as documented by our analysis in comparison to previously reported studies.
3. Detection of vsRNAs by deep sequencing indicates the activation of the RNAi response to both persistent and pathogenic infection of BmCPV.  S3 Dataset. Expression data of the 308 genes that are differentially expressed following pathogenic infection in both library pairs obtained by deep sequencing (2c/2inf and 4c/4inf pairs, corresponding to persistent/pathogenic infection at 2 nd and 4 th instar stages). Shown are RPKMs for each gene in each individual library (after trimming), as well as their average RPKMs in persistently and pathogenically infected samples. Total (not normalized) reads for each gene are also shown. The last two columns show the differential expression in pathogenically versus persistently infected midgut tissue for 2 nd and 4 th instar developmental stages. Criteria for the selection of genes are described in detail in Materials and Methods section. The graphs depict mean values of expression normalized to the housekeeping gene actin 3, as measured for two biological and two technical replicates (+SE). For clarity, four different graphs (a-d) of relative mRNA levels are shown, in which genes with similar mRNA levels are grouped (see different scales in the graphs). High error bars obtained for expression of genes may reflect differences between samples in exact developmental stage or progression of viral infection. Expression of genes may be more strictly developmentally regulated or be more sensitive to progression of viral infection. See Table 2  Expression of several RNAi-related genes was analyzed by qRT-PCR in midgut samples of persistently and pathogenically infected 2 nd instar larvae. The graphs depict mean values of expression normalized to the housekeeping gene actin 3, as measured for two biological and two technical replicates (+SE). For clarity, four different graphs (a-d) of relative mRNA levels are shown, in which genes with similar mRNA levels are grouped (see different scales in the graphs). High error bars obtained for expression of genes may reflect differences between samples in exact developmental stage or progression of viral infection. Expression of genes may be more strictly developmentally regulated or more sensitive to progression of viral infection. See Table 6 and [64] for further explanation on gene identity and function. (TIF) S1 Table. Sequences of primers used in qRT-PCR reactions for mRNA expression quantification. Genes are identified by their SilkDB ID (www.silkdb.org) and their GenBank Accession Number (www.ncbi.nlm.nih.gov/). (DOCX)