Journal of Invertebrate Pathology 115 (2014) 8–13
Contents lists available at ScienceDirect
Journal of Invertebrate Pathology journal homepage: www.elsevier.com/locate/jip
Digital Gene Expression analysis in the midgut of 4008 silkworm strain infected with cytoplasmic polyhedrosis virus Kun Gao a,b,c, Xiang-yuan Deng a, He-ying Qian b,c, Guangxing Qin b,c, Xi-jie Guo b,c,⇑ a
College of Biotechnology and Chemical Engineering, Jiangsu University of Science and Technology, Zhenjiang 212003, Jiangsu, China Sericultural Research Institute, Chinese Academy of Agricultural Sciences, Zhenjiang 212018, Jiangsu, China c Sericultural Research Institute, Jiangsu University of Science and Technology, Zhenjiang 212018, Jiangsu, China b
a r t i c l e
i n f o
Article history: Received 16 September 2013 Accepted 28 October 2013 Available online 6 November 2013 Keywords: Bombyx mori Digital Gene Expression Cytoplasmic polyhedrosis virus Midgut
a b s t r a c t Digital Gene Expression was performed to investigate the midgut transcriptome profile of 4008 silkworm strain orally infected with BmCPV. A total of 4,498,263 and 4,258,240 clean tags were obtained from the control and BmCPV-infected larvae. A total of 752 differentially expressed genes were detected, of which 649 were upregulated and 103 were downregulated. Analysis results of the Kyoto Encyclopedia of Genes and Genomes pathway showed that 334 genes were involved in the ribosome and RNA transport pathways. Moreover, 408 of the 752 differentially expressed genes have a GO category and can be categorized into 41 functional groups according to molecular function, cellular component and biological process. Differentially expressed genes involved in signaling, gene expression, metabolic process, cell death, binding, and catalytic activity changes were detected in the expression profiles. Quantitative real-time PCR was performed to verify the expression of these genes. The upregulated expression levels of Calreticulin, FK506-binding protein, and protein kinase c inhibitor gene probably led to a calcium-dependent apoptosis in the BmCPV-infected cells. The results of this study may serve as a basis for future research not only on the molecular mechanism of BmCPV invasion but also on the anti-BmCPV mechanism of silkworm. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Bombyx mori cytoplasmic polyhedrosis virus (BmCPV), which belongs to the CPV subfamily within the genus Cypovirus of family Reoviridae, is an important viral pathogen in silkworms (Hill et al., 1999). BmCPV specifically infects epithelial cells in the midgut of silkworms, causing white wrinkles in the posterior part of the midgut, a typical symptom of the disease (Magnoler, 1974). Studies on insect antiviral system and interactive mechanisms between the host cells and BmCPV are still at their infancy and remain deficient. A microarray analysis was carried out to examine the gene transcript changes in the midgut of BmCPV-infected and normal p50 silkworm larvae (Wu et al., 2011). At 72 h post-inoculation, 258 genes exhibited at least 2.0-fold differences in expression levels. The expression levels of genes related to the degradation of valine, leucine, and isoleucine, as well as to retinol metabolism and vitamin B6 metabolism, were downregulated. Meanwhile, the expression levels of genes involved in ribosomal and proteasomal pathways were upregulated. These results suggest that BmCPV
⇑ Corresponding author at: Sericultural Research Institute, Jiangsu University of Science and Technology, Zhenjiang 212018, Jiangsu, China. Fax: +86 511 85628183. E-mail address: [email protected] (X.-j. Guo). 0022-2011/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jip.2013.10.016
infection disrupts protein and amino acid metabolism and induces a series of major physiologic and pathologic changes in silkworms. Recent advances in the development of ultra high-throughput deep sequencing technologies are making a huge impact on genomic research. These next generation sequencing systems, such as the Solexa/Illumina Genome Analyzer and the ABI/SOLiD Gene Sequencer, can sequence in parallel millions of DNA molecules derived directly from mRNA without using bacterial clones (Cloonan and Grimmond, 2008; Morozova and Marra, 2008; Wang et al., 2009). The direct sequencing yields libraries of short sequences (25–50 nucleotides), referred to as RNA-Seq data or Digital Gene Expression (DGE) data. Sequencing-based methods generate absolute rather than relative gene expression measurements and avoid many of the inherent limitations of microarray analysis which has been the most commonly used technology for transcriptome profiling over the last decade. In order to clarify the molecular mechanism of BmCPV infection to silkworm and host response, we examine the transcriptome profiles of a susceptible silkworm strain 4008 to screen and detect more disease-related differentially expressed genes by employing the second generation Illumina Genome Analyzer platform (GA II) to perform a SAGE-derived Digital Gene Expression (DGE) analysis. Furthermore, some selected differentially expressed genes were verified by qRT-PCR.
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13
2. Materials and methods 2.1. Silkworm rearing and Virus inoculation The silkworm strain 4008 used in this study was provided by the National Center for Silkworm Genetic Resources Preservation of the Chinese Academy of Agricultural Sciences. The larvae were reared at standard temperature and under a photoperiod of 12 h of light and 12 h of dark up to second molting for virus inoculation. BmCPV viral stock was suspended in distilled water at a concentration of 108 polyhedra per ml. 1 ll virus suspension was fed to newly molted third instar larvae of silkworm. The infection dose was calculated to be 1 105 polyhedra per larva. The control larvae were fed with 1 ll distilled water. 2.2. DGE tag profiling We analyzed the gene expression of the 4008 at 72 h in the fifth instar after orally infected by BmCPV using the DGE method. Considering the effects of individual differences, we dissected 15 larvae per group as a sample for DGE, respectively. The mRNA was purified from 6 lg of total RNA from each sample, and then used for synthesizing the first- and second-strand cDNA. NlaIII, a 4 base-recognition enzyme, was used to digest this cDNA, and Illumina adaptor 1 was ligated. Mmel was used to digest at a region 17 bp downstream of the CATG site; afterwards, an Illumina adaptor 2 was ligated at the 30 end. After 15 cycles of linear PCR amplification, 105-bp fragments were purified by 6% TBE PAGE gel electrophoresis. After denaturation, the single-chain molecules were fixed onto an Illumina Sequencing Chip (flowcell). Each molecule was grown into a single-molecule cluster sequencing template through in situ amplification. Subsequently, four types of nucleotides, which were labeled by one of four colors, were added to the chip, and then sequencing was performed with the sequencing by synthesis method. Each tunnel could generate millions of raw reads, each with a sequencing length of 49 bp. Adaptor sequences, low quality tags (tags with unknown nucleotides N), empty reads (reads with only 30 adaptor sequences, but no tags), tags that were too short or too long, and tags with only one copy (probable sequencing error) were filtered to obtain clean tags. A preprocessed database of all possible CATG + 17-nucleotide tag sequences was created using the transcriptome reference database generated from SilkDB v2.0 database. For annotation, all clean tags were mapped to the reference sequence, and a mismatch of only 1 bp was considered. Clean tags mapped to the reference sequences from multiple genes were filtered. The remaining clean tags were designated as unambiguous clean tags. For gene expression analysis, the number of unambiguous clean tags for each gene was calculated, and then normalized to the number of transcripts per million (TPM) clean tags. 2.3. Analysis of differential gene expression The differential expression detection of genes across samples was performed using a rigorous algorithm method (Benjamini and Yekutieli, 2001). False discovery rate (FDR) was used to determine the P value threshold in multiple tests and analyses. We obtained the significance of the gene expression difference through an FDR 60.001 and the absolute value of log 2 ratio P1. 2.4. GO and KEGG The gene ontology (GO) classification system was used to determine the possible functions of all differentially expressed genes. P value was calculated by GO (http://www.geneontology.org/)
9
and Bonferoni corrected. A corrected P value 60.05 was selected as a threshold for significant enrichment of the gene sets. WEGO (Web Gene Ontology Annotation Plot) software is used for visualizing, comparing and plotting GO annotation results (Ye et al., 2006). Pathway enrichment analysis can further identify significantly enriched metabolic pathways or signal transduction pathways using the KEGG database. Pathways with a Q value 60.05 are designated as significantly enriched pathways in DGEs. 2.5. qRT-PCR analysis The genes selected according to the DGE-tag copy number were evaluated, and some of these genes were investigated by qRT-PCR. Bmactin A3 was used as a reference gene and the primers used in this study were in Table S1. 1 lg of total RNA from each sample was used to synthesize the first strand cDNA using the PrimeScript Reverse Transcriptase kit (TaKaRa) according to the protocol of the manufacturer. qRT-PCR was carried out in an ABI PRISMÒ 7300 Sequence Detection System (Applied Biosystems) using SYBR Green Supermix (TaKaRa) according the instructions of the manufacturer. The thermal cycle conditions were 94 °C for 10 min for denaturation, followed by 40 cycles of 94 °C for 15 s, and 60 °C for 31 s for annealing, and then an extension. The expression level of genes was calculated by 2DDCt method (Livad and Schmittgen, 2001), and the value stood for an nfold difference relative to the calibrator. All data were given in terms of relative mRNA expression as mean ± SE. 3. Results 3.1. Confirmation of infection At about 72 h post-inoculation of BmCPV, all the inoculated silkworm larvae were infected and diseased, exhibiting white wrinkles on the midgut as a typical symptom. Further confirmation was performed using a microscope by observing the BmCPV polyhedra. 3.2. Analysis of DGE libraries Two DGE libraries were made from midgut of the control and BmCPV-infected 4008 silkworm using Solexa (Illumina) technology. Total 4,660,291 and 4,389,132 raw tags were generated in the control and BmCPV infected DGE libraries, respectively. A total of 4,498,263 and 4,258,240 clean tags corresponding to 49,365 and 60,834 distinct clean tags were filtered from the raw tags (Table 1). The distribution of the total and distinct tags over the different tag abundance categories showed highly similar patterns in each DGE library. More than 85% of the total clean tags had a copy number higher than 100 counts, whereas less than 7% of distinct tags had a copy number higher than 100 counts (Fig. 1). The sequencing saturation was analyzed to estimate whether the sequencing depth was sufficient for transcriptome coverage. The results showed that the number of detected genes was almost saturated when the total tag number reached 2 million or higher. Our sequencing depths reached approximately 3.5 million in each DGE library (Fig. S1), which satisfied the requirement for the experiment. So the two DGE libraries were reliable. 3.3. Analysis of tag mapping We mapped the tag sequences of the two DGE libraries to the reference database of silkworm, which contains 16,329 reference genes including 13,328 genes with the CATG site. In the control and BmCPV-infected DGE library, 31.92% and 29.61% of the clean tags corresponding to 20.3% and 19.71% of distinct tags were
10
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13
Table 1 Distributions of tags in two DGE libraries. Category
Parameter
Control (4008)
BmCPV-infected (4008)
Raw data
Total Distinct tags Total number Distinct tag number Total number Total% of clean tags Distinct tag number Distinct Tag% of clean Total number Total% of clean tags Distinct tag number Distinct Tag% of clean Number % of ref genes Number % of ref genes Total number Total% of clean tags Distinct tag number Distinct Tag% of clean Total number Total% of clean tags Distinct tag number Distinct Tag% of clean
4660291 122861 4498263 49365 1435964 31.92% 10023 20.3% 792073 17.61% 9387 19.02% 3956 24.23% 3728 22.83% 1547514 34.4% 17254 34.95% 840616 18.69% 14132 28.63%
4389132 144664 4258240 60834 1260841 29.61% 11991 19.71% 851711 20.00% 11372 18.69% 4621 28.3% 4359 26.69% 1739966 40.86% 22582 37.12% 755891 17.75% 16145 26.54%
Clean tags
All tags mapping to gene
Unambiguous tags Mapping to gene All tag-mapped genes Unambiguous Tag-mapped genes Mapping to genome
Unknown tags
tags
tags
tags
tags
Fig. 1. Distribution of total clean tags and distinct clean tags in each DGE library. C3: the control midgut of 4008 silkworm of the 3rd day in the fifth-instar; T3: the BmCPVinfected midgut of 4008 silkworm of the 3rd day in the fifth-instar. The numbers in square brackets indicate the range of copy numbers of each tag category. The data in parentheses indicate the number and percentage of corresponding tags among the total clean tags and distinct clean tags.
mapped to a gene in the reference database, 34.4% and 40.86% of the clean tags could be mapped to genome of silkworm, while 18.69% and 17.75% of the clean tags were unknown tags (Table 1). The number of unambiguous tags for each gene was normalized to TPM to accurately measure the gene expression level. The results showed that the mRNA transcribed from the major types of genes were represented in low copy number, and only a small proportion of genes were highly expressed (Fig. S2). 3.4. GO analysis and KEGG pathways A total of 752 differentially expressed genes were detected, of which 648 were upregulated and 103 were downregulated
(Fig. 2). To understand the functions of these differentially expressed genes, all of the 752 differentially expressed genes were mapped to terms in GO database and compared with the whole genome background. GO has three ontologies, namely, molecular function, cellular component, and biological process. 408 differentially expressed genes including 354 up-regulated genes and 54 down-regulated genes have a GO ID and can be categorized into a total of 41 functional groups using the WEGO software (Fig. 3). In each of the three main categories of the GO classification, ‘‘cellular process’’, ‘‘metabolic process’’, ‘‘cell and cell parts’’, ‘‘catalytic activity’’, ‘‘binding’’ terms are dominant. Interestingly, some genes related to the cell-proliferation, death, growth, locomotion, signaling, membrane-enclosed lumen, structural molecule, biological
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13
11
4986 significantly changed tags were found and mapped to 752 genes with 648 up-regulated and 103 down-regulated genes in DGE libraries. A total of 21 differentially expressed genes including 11 up-regulated and 10 down-regulated genes from DGE libraries were selected for qRT-PCR analysis to validate the DGE data. The results showed that 20 genes were demonstrated to have a concordant direction of change for both DGE and qRT-PCR except Bm_nscaf463_32 gene (Fig. 4).
4. Discussion
Fig. 2. Differentially expressed genes. The red part represents those genes upregulated in BmCPV-infected group compared to control group. The green part shows those genes down-regulated in BmCPV-infected group. The blue part shows those genes without expression difference between these two samples. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
process, nucleic acid binding, gene expression were all up-regulated in BmCPV-infected midgut of 4008 silkworm. KEGG (http://www.genome.jp/kegg) ontology assignments were used to classify the functional annotations of the identified genes to further understand the biological functions of the genes. Among the differentially expressed genes, 613 were mapped to 212 pathways in the KEGG database in 4008 DGE class. 19 terms was significantly enriched (Q value 60.05) (Table S2). 3.5. Identification and verification of differentially expressed genes The differentially expressed tags between the two samples were identified by FDR 60.001 and |log 2 ratio|P1. A total of
Bombyx mori cytplasmic polyhedrosis virus (BmCPV) is an important pathogen to the silkworm, always causing serious damage to sericulture production. However, the molecular mechanism of infection and defense response of the silkworm remains unclear. Silkworm strains p50 and 4008 are different in resistance to BmCPV infection according to previous studies. It showed that p50 is more resistant and 4008 is more susceptible (Xu et al., 2002; Wu et al., 2013). The previous microarray study by Wu et al. (2011) was focused on strain p50 but our present study was focused on the susceptible strain 4008. Two DGE libraries were made from midgut of the control and BmCPV-infected 4008 silkworm using Solexa (Illumina) technology with 4,660,291 and 4,389,132 raw tags, respectively. A total of 752 differentially expressed genes were detected, of which 648 were upregulated and 103 were downregulated (Fig. 2). The number of upregulated genes was markedly higher (over 6-fold) than the number of downregulated genes, which is different to the previous microarray study of B. mori larvae p50 strain infected with BmCPV with 135 upregulated and 123 downregulated genes in hybridization-based (microarrays) technology. In fact, we also examined the transcriptome profiling between the infected and mock-infected larvae of strain p50 with the same Digital Gene Expression approach and detected not only less number of differentially expressed genes but also less number of up-regulated genes than down-regulated ones, which is consistent to the result of the previous microarray study (data not provided). Therefore, the differential expression of genes in the two silkworm strains might be
Fig. 3. Histogram presentation of Gene Ontology classification. The results are summarized in three main categories: biological process, cellular component and molecular function. The y-axis on the left indicates the percentage of a specific category of genes. The y-axis on the right indicates the number of genes in a category.
12
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13 15 gene expression in qRT-PCR gene expression in DGE
Relative expression level (Fold change=log2 Ratio)
10 5
0
-5 -10
-15
Genes Fig. 4. Differential expression levels of candidate genes in BmCPV-infected and control midgut of 4008 silkworm larvae. The y-axis indicates the relative expression level of candidate gene mRNA transcripts (Fold change = log 2 Ratio). The x-axis indicates the candidate genes (from left to right) with the names of Bm_nscaf1690_036, Bm_nscaf2829_122, Bm_nscaf2964_026, Bm_nscaf2902_170, Bm_nscaf3031_071, Bm_nscaf2891_038, Bm_nscaf463_32, Bm_nscaf3062_05, Bm_nscaf2204_100, Bm_nscaf2529_090, Bm_nscaf2204_124, Bm_nscaf2829_149, Bm_nscaf2828_137, Bm_nscaf3075_32, Bm_nscaf98_23, Bm_nscaf2901_47, Bm_nscaf3099_033, Bm_nscaf3058_068, Bm_nscaf2575_003, Bm_nscaf2887_005 and Bm_nscaf2529_052. The genes with different expression regulation by qRT-PCR were indicated with diamonds ( ). Vertical bars represent the mean ± SE (n = 3), and bars in different time marked with star were significantly different from each other (p < 0.01).
involved in different resistance or defense response to BmCPV infection, which need to be clarified in our future research. Among the differentially expressed genes, 613 genes were mapped to 212 pathways in the KEGG database. It showed that the ribosome and RNA transport pathways were significantly enriched. Genes related to the terms ‘‘metabolic process,’’ ‘‘structural molecule activity,’’ and ‘‘cell or cell part’’ were dominant in the BmCPV-infected 4008 silkworm. Meanwhile, genes related to ribosome, gene expression, translation, and ribonucleoprotein complex biogenesis were all upregulated (Fig. 3). This result is mostly similar with the result of a DNA microarray employed by Wu et al. (2011). It suggested that BmCPV infection extensively affects known gene expression pathways. During the period of BmCPV infection, a large number of genes used for translation were overexpressed to meet the needs of the enormous proliferation of BmCPV. Among the 752 differentially expressed genes, 11 upregulated and 10 downregulated genes were verified via qRT-PCR (Table S1). Among the 20 concordant genes, only 12 genes, Bm_nscaf1690_036, Bm_nscaf2964_026, Bm_nscaf3031_071, Bm_nscaf3062_05, Bm_nscaf2204_100, Bm_nscaf2529_090, Bm_nscaf2829_149, Bm_nscaf98_23, Bm_nscaf2901_47, Bm_nscaf3058_068, Bm_nscaf2575_003, Bm_nscaf2887_005 were identified with FDR 60.001 and |log 2 ratio|P1 verified using qRT-PCR (Fig. 4). Lepidopteron species can resist baculovirus infection through selective apoptosis of infected cells in the midgut epithelium and through sloughing off infected cells (Clarke and Clem, 2002). Thus, the upregulated expression of the programmed cell death (PCD) protein 5 gene (Bm_nscaf3031_071) in the midgut cells may induce apoptosis in BmCPV-infected silkworm (Chen et al., 2013, 2001; Wang et al., 2007). Meanwhile, PCD is mediated via a nongenomic pathway that includes Ca(2+) as a second messenger and the activation of protein kinase C/caspase-3-like protease (Manaboon et al., 2009); Bm_nscaf1681_279, which encodes a Calreticulin (CRT) gene, is one of the major Ca(2+) binding chaperone proteins of the endoplasmic reticulum (ER) and plays a critical role in the activation of Ca(2+)-dependent apoptosis (Lim et al., 2008). In the cells over-expressing CRT, the intracellular calcium concentration was significantly increased and the activity of protein kinase c mRNA were reduced. It suggested that the up-regulated expression of CRT can disrupts intracellular calcium regulation leading to a calcium-dependent apoptosis in the BmCPV-infected cells. Furthermore, Cytosolic Ca(2+) concentration is crucially affected by the activity of the intracellular store (the sarcoplasmic reticulum, SR), which regulates Ca(2+) release (Berridge et al., 2003; McCarron et al., 2004). Two distinct classes of Ca(2+) release
channels were the inositol 1,4,5-trisphosphate (IP3) receptor and the ryanodine receptor. The important modulatory proteins that interact with these channels were FK506 binding proteins (FKBPs), including FKBP12 and its isoform-FKBP12.6 (Van Acker et al., 2004; Zissimopoulos et al., 2007). These regulatory proteins are reportedly important modulators of intracellular Ca(2+) release. In addition, Cyclic ADP-ribose (cADPR) binds to FKBP 12.6 on rat islet ryanodine receptor and that the binding of cADPR to FKBP12.6 frees the ryanodine receptor from FKBP12.6, causing it to release Ca(2+) (Noguchi et al., 1997). In addition, recent study have showed that the Ca(2+) entry channels in prothoracic gland cells of B. mori are probably modulated by protein kinases A and C. Application of protein kinases A and C inhibitors to the prothoracic glands cells of the silkworm, B. mori, resulted in slow and gradual increases in intracellular Ca(2+) (Dedos et al., 2008). FK506-binding protein (FKBP), which has also been identified as an inhibitor of protein kinase C (Nakazawa et al., 2005) may play a role in slow and gradual increases in intracellular Ca(2+). So the up-regulated expression of FK506-binding protein gene (Bm_nscaf1898_230) and protein kinase c inhibitor gene (Bm_nscaf3090_7) probably induced the increase of intracellular Ca(2+) in the BmCPV-infected cells and triggered the cell death. Early and pivotal events in apoptosis are now known to occur in the mitochondria and the endoplasmic reticulum(ER), and release of cytochrome c from the mitochondria and calcium from the endoplasmic recticulum into the cytosol is considered to be requisites for apoptosis in response to different stimuli (Wang and ElDeiry, 2004). FK506-binding protein, which has been identified as an inhibitor of protein kinase C is also an important modulatory protein that interacts with Ca(2+) release channels and may play a role in slow and gradual increases in intracellular Ca(2+). Calreticulin (CRT) is one of the major Ca2+ binding chaperone proteins of the ER and an unusual luminal ER protein. In the cells overexpressing CRT, the intracellular calcium concentration was significantly increased and the activity of PKC was reduced, leading to calcium-dependent apoptosis. The function of these two genes in the calcium-dependent apoptosis in the silkworm and furthermore the host response to BmCPV by the calcium-dependent apoptosis should be clarified in our future study. Genes related to activity of digestive enzyme were all down-regulated in the BmCPV-infected silkworms. Those genes include, for example, a chymotrypsin-like serine protease (Bm_nscaf2901_47), an acidic lipase (Bm_nscaf2829_149) and an lipase-1 (Bm_nscaf2529_052) (Fig. 4) etc. BmCPV is specific for epithelial cells in the silkworm midgut, which acts as the main site for virus
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13
multiplication (Miyajima and Kawase, 1968). As disease progresses, white wrinkles typically occur in the posterior part of the midgut, which may affect the digestive and the absorptive functions of the midgut. So the fatty acid binding protein gene (Bm_nscaf3058_068) and sugar transporter gene (Bm_nscaf2575_003) were also downregulated in the BmCPV-infected silkworm. In conclusion, we obtained 752 differentially expressed genes from the control and BmCPV-infected midgut of 4008 silkworm through DGE analysis. Analysis results of the KEGG pathway showed that the ribosome and RNA transport pathways were significantly enriched. The expression profiles of some important differentially expressed genes involved in signaling, gene expression, metabolic process, cell death, binding, and catalytic activity changed. However, most of these differentially expressed genes were up-regulated. These results suggest that the upregulated expression of CRT, FK506-binding protein, and protein kinase c inhibitor gene probably led to a calcium-dependent apoptosis in the BmCPV-infected cells. Therefore, the death or apoptosis program may be an antiviral response to defend the host from BmCPV proliferation. But further studies with other techniques are still needed, Such as over-expression, RNAi, western blotting, co-immunoprecipitation, and so on. The current results provide important clues for future research on the molecular mechanism of BmCPV infection. Acknowledgments This work was supported by the Jiangsu Postdoctoral Science Foundation of China (No. 1202026C) and the National Basic Research Program of China (Grant No. 2012CB114600). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jip.2013.10.016. References Benjamini, Y., Yekutieli, D., 2001. The control of the false discovery rate in multiple 348 testing under dependency. Ann. Stat. 29, 1165–1188. Berridge, M.J., Bootman, M.D., Roderick, H.L., 2003. Calcium signalling: dynamics, homeostasis and remodelling. Nat. Rev. Mol. Cell Biol. 4, 517–529. Chen, Y., Sun, R., Han, W., Zhang, Y., Song, Q., Di, C., Ma, D., 2001. Nuclear translocation of PDCD5 (TFAR19): an early signal for apoptosis? FEBS Lett. 509, 191–196. Chen, C.H., Jiang, Z., Yan, J.H., Yang, L., Wang, K., Chen, Y.Y., Han, J.Y., Zhang, J.H., Zhou, C.M., 2013. The Involvement of Programmed Cell Death 5 (PDCD5) in the regulation of apoptosis in cerebral ischemia/reperfusion injury. CNS Neurosci. Ther. 19, 566–576. Clarke, T.E., Clem, R.J., 2002. Lack of involvement of haemocytes in the establishment and spread of infection in Spodoptera frugiperda larvae
13
infected with the baculovirus Autographa californica M nucleopolyhedrovirus by intrahaemocoelic injection. J. Gen. Virol. 83, 1565–1572. Cloonan, N., Grimmond, S.M., 2008. Transcriptome content and dynamics at singlenucleotide resolution. Genome Biol. 9, 234. Dedos, S.G., Kaltofen, S., Birkenbeil, H., 2008. Protein kinase A and C are ‘‘Gatekeepers’’ of capacitative Ca2+ entry in the prothoracic gland cells of the silkworm, Bombyx mori. J. Insect Physiol. 54, 878–882. Hill, C.L., Booth, T.F., Stuart, D.I., Mertens, P.P., 1999. Lipofectin increases the specific activity of cypovirus particles for cultured insect cells. J. Virol. Methods 78, 177– 189. Livad, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25, 402–408. Lim, S., Chang, W., Lee, B.K., Song, H., Hong, J.H., Lee, S., Song, B.W., Kim, H.J., Cha, M.J., Jang, Y., Chung, N., Choi, S.Y., Hwang, K.C., 2008. Enhanced calreticulin expression promotes calcium-dependent apoptosis in postnatal cardiomyocytes. Mol. Cells 25, 390–396. Magnoler, A., 1974. Effects of a cytoplasmic polyhedrosis on larval and postlarval stages of the gypsy moth, Porthetria dispar. J. Invertebr. Pathol. 23, 263–274. Manaboon, M., Iga, M., Iwami, M., Sakurai, S., 2009. Intracellular mobilization of Ca2+ by the insect steroid hormone 20-hydroxyecdysone during programmed cell death in silkworm anterior silk glands. J. Insect Physiol. 55, 122–128. McCarron, J.G., Bradley, K.N., MacMillan, D., Chalmers, S., Muir, T.C., 2004. The sarcoplasmic reticulum, Ca2+ trapping, and wave mechanisms in smooth muscle. News Physiol. Sci. 19, 138–147. Miyajima, S., Kawase, S., 1968. Changes in virus-infectivity titer in the hemolymph and midgut during the course of a cytoplasmic polyhedrosis in the silkworm. J. Invertebr. Pathol. 12, 329–334. Morozova, O., Marra, M.A., 2008. Applications of next-generation sequencing technologies in functional genomics. Genomics 92, 255–264. Nakazawa, T., Takasawa, S., Noguchi, N., Nata, K., Tohgo, A., Mori, M., Nakagawara, K., Akiyama, T., Ikeda, T., Yamauchi, A., Takahashi, I., Yoshikawa, T., Okamoto, H., 2005. Genomic organization, chromosomal localization, and promoter of human gene for FK506-binding protein 12.6. Gene 360, 55–64. Noguchi, N., Takasawa, S., Nata, K., Tohgo, A., Kato, I., Ikehata, F., Yonekura, H., Okamoto, H., 1997. Cyclic ADP-ribose binds to FK506-binding protein 12.6 to release Ca2+ from islet microsomes. J. Biol. Chem. 272, 3133–3136. Van Acker, K., Bultynck, G., Rossi, D., Sorrentino, V., Boens, N., Missiaen, L., De Smedt, H., Parys, J.B., Callewaert, G., 2004. The 12 kDa FK506-binding protein, FKBP12, modulates the Ca(2+)-flux properties of the type-3 ryanodine receptor. J. Cell Sci. 117, 1129–1137. Wang, S., El-Deiry, W.S., 2004. Cytochrome c: a crosslink between the mitochondria and the endoplasmic reticulum in calcium-dependent apoptosis. Cancer Biol. Ther. 3, 44–46. Wang, N., Lu, H.S., Guan, Z.P., Sun, T.Z., Chen, Y.Y., Ruan, G.R., Chen, Z.K., Jiang, J., Bai, C.J., 2007. Involvement of PDCD5 in the regulation of apoptosis in fibroblast-like synoviocytes of rheumatoid arthritis. Apoptosis 12, 1433–1441. Wang, Z., Gerstein, M., Snyder, M., 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57–63. Wu, P., Qian, H.Y., Xu, A.Y., Chen, T., Hou, Q.R., Guo, X.J., Li, L., 2013. Analysis on resistance of various Bombyx mori Varieties to cytoplasmic polyhedrosis virus infection by quantitative PCR. Sci. Sericul. 39, 0289–0294. Wu, P., Wang, X., Qin, G.X., Liu, T., Jiang, Y.F., Li, M.W., Guo, X.J., 2011. Microarray analysis of the gene expression profile in the midgut of silkworm infected with cytoplasmic polyhedrosis virus. Mol. Biol. Rep. 38, 333–341. Xu, A.Y., Li, M.W., Lin, C.Q., Zhang, Y.H., Hou, C.X., 2002. The comparision of the resistibility to Bombyx mori cytoplasmic polyhedrosis virus among silkworm germplasm resourse. Sci. Sericul. 28, 157–159. Ye, J., Fang, L., Zheng, H., Zhang, Y., Chen, J., Zhang, Z., Wang, J., Li, S., Li, R., Bolund, L., 2006. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 34, W293–W297. Zissimopoulos, S., Docrat, N., Lai, F.A., 2007. Redox sensitivity of the ryanodine receptor interaction with FK506-binding protein. J. Biol. Chem. 282, 6976–6983.
Contents lists available at ScienceDirect
Journal of Invertebrate Pathology journal homepage: www.elsevier.com/locate/jip
Digital Gene Expression analysis in the midgut of 4008 silkworm strain infected with cytoplasmic polyhedrosis virus Kun Gao a,b,c, Xiang-yuan Deng a, He-ying Qian b,c, Guangxing Qin b,c, Xi-jie Guo b,c,⇑ a
College of Biotechnology and Chemical Engineering, Jiangsu University of Science and Technology, Zhenjiang 212003, Jiangsu, China Sericultural Research Institute, Chinese Academy of Agricultural Sciences, Zhenjiang 212018, Jiangsu, China c Sericultural Research Institute, Jiangsu University of Science and Technology, Zhenjiang 212018, Jiangsu, China b
a r t i c l e
i n f o
Article history: Received 16 September 2013 Accepted 28 October 2013 Available online 6 November 2013 Keywords: Bombyx mori Digital Gene Expression Cytoplasmic polyhedrosis virus Midgut
a b s t r a c t Digital Gene Expression was performed to investigate the midgut transcriptome profile of 4008 silkworm strain orally infected with BmCPV. A total of 4,498,263 and 4,258,240 clean tags were obtained from the control and BmCPV-infected larvae. A total of 752 differentially expressed genes were detected, of which 649 were upregulated and 103 were downregulated. Analysis results of the Kyoto Encyclopedia of Genes and Genomes pathway showed that 334 genes were involved in the ribosome and RNA transport pathways. Moreover, 408 of the 752 differentially expressed genes have a GO category and can be categorized into 41 functional groups according to molecular function, cellular component and biological process. Differentially expressed genes involved in signaling, gene expression, metabolic process, cell death, binding, and catalytic activity changes were detected in the expression profiles. Quantitative real-time PCR was performed to verify the expression of these genes. The upregulated expression levels of Calreticulin, FK506-binding protein, and protein kinase c inhibitor gene probably led to a calcium-dependent apoptosis in the BmCPV-infected cells. The results of this study may serve as a basis for future research not only on the molecular mechanism of BmCPV invasion but also on the anti-BmCPV mechanism of silkworm. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Bombyx mori cytoplasmic polyhedrosis virus (BmCPV), which belongs to the CPV subfamily within the genus Cypovirus of family Reoviridae, is an important viral pathogen in silkworms (Hill et al., 1999). BmCPV specifically infects epithelial cells in the midgut of silkworms, causing white wrinkles in the posterior part of the midgut, a typical symptom of the disease (Magnoler, 1974). Studies on insect antiviral system and interactive mechanisms between the host cells and BmCPV are still at their infancy and remain deficient. A microarray analysis was carried out to examine the gene transcript changes in the midgut of BmCPV-infected and normal p50 silkworm larvae (Wu et al., 2011). At 72 h post-inoculation, 258 genes exhibited at least 2.0-fold differences in expression levels. The expression levels of genes related to the degradation of valine, leucine, and isoleucine, as well as to retinol metabolism and vitamin B6 metabolism, were downregulated. Meanwhile, the expression levels of genes involved in ribosomal and proteasomal pathways were upregulated. These results suggest that BmCPV
⇑ Corresponding author at: Sericultural Research Institute, Jiangsu University of Science and Technology, Zhenjiang 212018, Jiangsu, China. Fax: +86 511 85628183. E-mail address: [email protected] (X.-j. Guo). 0022-2011/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jip.2013.10.016
infection disrupts protein and amino acid metabolism and induces a series of major physiologic and pathologic changes in silkworms. Recent advances in the development of ultra high-throughput deep sequencing technologies are making a huge impact on genomic research. These next generation sequencing systems, such as the Solexa/Illumina Genome Analyzer and the ABI/SOLiD Gene Sequencer, can sequence in parallel millions of DNA molecules derived directly from mRNA without using bacterial clones (Cloonan and Grimmond, 2008; Morozova and Marra, 2008; Wang et al., 2009). The direct sequencing yields libraries of short sequences (25–50 nucleotides), referred to as RNA-Seq data or Digital Gene Expression (DGE) data. Sequencing-based methods generate absolute rather than relative gene expression measurements and avoid many of the inherent limitations of microarray analysis which has been the most commonly used technology for transcriptome profiling over the last decade. In order to clarify the molecular mechanism of BmCPV infection to silkworm and host response, we examine the transcriptome profiles of a susceptible silkworm strain 4008 to screen and detect more disease-related differentially expressed genes by employing the second generation Illumina Genome Analyzer platform (GA II) to perform a SAGE-derived Digital Gene Expression (DGE) analysis. Furthermore, some selected differentially expressed genes were verified by qRT-PCR.
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13
2. Materials and methods 2.1. Silkworm rearing and Virus inoculation The silkworm strain 4008 used in this study was provided by the National Center for Silkworm Genetic Resources Preservation of the Chinese Academy of Agricultural Sciences. The larvae were reared at standard temperature and under a photoperiod of 12 h of light and 12 h of dark up to second molting for virus inoculation. BmCPV viral stock was suspended in distilled water at a concentration of 108 polyhedra per ml. 1 ll virus suspension was fed to newly molted third instar larvae of silkworm. The infection dose was calculated to be 1 105 polyhedra per larva. The control larvae were fed with 1 ll distilled water. 2.2. DGE tag profiling We analyzed the gene expression of the 4008 at 72 h in the fifth instar after orally infected by BmCPV using the DGE method. Considering the effects of individual differences, we dissected 15 larvae per group as a sample for DGE, respectively. The mRNA was purified from 6 lg of total RNA from each sample, and then used for synthesizing the first- and second-strand cDNA. NlaIII, a 4 base-recognition enzyme, was used to digest this cDNA, and Illumina adaptor 1 was ligated. Mmel was used to digest at a region 17 bp downstream of the CATG site; afterwards, an Illumina adaptor 2 was ligated at the 30 end. After 15 cycles of linear PCR amplification, 105-bp fragments were purified by 6% TBE PAGE gel electrophoresis. After denaturation, the single-chain molecules were fixed onto an Illumina Sequencing Chip (flowcell). Each molecule was grown into a single-molecule cluster sequencing template through in situ amplification. Subsequently, four types of nucleotides, which were labeled by one of four colors, were added to the chip, and then sequencing was performed with the sequencing by synthesis method. Each tunnel could generate millions of raw reads, each with a sequencing length of 49 bp. Adaptor sequences, low quality tags (tags with unknown nucleotides N), empty reads (reads with only 30 adaptor sequences, but no tags), tags that were too short or too long, and tags with only one copy (probable sequencing error) were filtered to obtain clean tags. A preprocessed database of all possible CATG + 17-nucleotide tag sequences was created using the transcriptome reference database generated from SilkDB v2.0 database. For annotation, all clean tags were mapped to the reference sequence, and a mismatch of only 1 bp was considered. Clean tags mapped to the reference sequences from multiple genes were filtered. The remaining clean tags were designated as unambiguous clean tags. For gene expression analysis, the number of unambiguous clean tags for each gene was calculated, and then normalized to the number of transcripts per million (TPM) clean tags. 2.3. Analysis of differential gene expression The differential expression detection of genes across samples was performed using a rigorous algorithm method (Benjamini and Yekutieli, 2001). False discovery rate (FDR) was used to determine the P value threshold in multiple tests and analyses. We obtained the significance of the gene expression difference through an FDR 60.001 and the absolute value of log 2 ratio P1. 2.4. GO and KEGG The gene ontology (GO) classification system was used to determine the possible functions of all differentially expressed genes. P value was calculated by GO (http://www.geneontology.org/)
9
and Bonferoni corrected. A corrected P value 60.05 was selected as a threshold for significant enrichment of the gene sets. WEGO (Web Gene Ontology Annotation Plot) software is used for visualizing, comparing and plotting GO annotation results (Ye et al., 2006). Pathway enrichment analysis can further identify significantly enriched metabolic pathways or signal transduction pathways using the KEGG database. Pathways with a Q value 60.05 are designated as significantly enriched pathways in DGEs. 2.5. qRT-PCR analysis The genes selected according to the DGE-tag copy number were evaluated, and some of these genes were investigated by qRT-PCR. Bmactin A3 was used as a reference gene and the primers used in this study were in Table S1. 1 lg of total RNA from each sample was used to synthesize the first strand cDNA using the PrimeScript Reverse Transcriptase kit (TaKaRa) according to the protocol of the manufacturer. qRT-PCR was carried out in an ABI PRISMÒ 7300 Sequence Detection System (Applied Biosystems) using SYBR Green Supermix (TaKaRa) according the instructions of the manufacturer. The thermal cycle conditions were 94 °C for 10 min for denaturation, followed by 40 cycles of 94 °C for 15 s, and 60 °C for 31 s for annealing, and then an extension. The expression level of genes was calculated by 2DDCt method (Livad and Schmittgen, 2001), and the value stood for an nfold difference relative to the calibrator. All data were given in terms of relative mRNA expression as mean ± SE. 3. Results 3.1. Confirmation of infection At about 72 h post-inoculation of BmCPV, all the inoculated silkworm larvae were infected and diseased, exhibiting white wrinkles on the midgut as a typical symptom. Further confirmation was performed using a microscope by observing the BmCPV polyhedra. 3.2. Analysis of DGE libraries Two DGE libraries were made from midgut of the control and BmCPV-infected 4008 silkworm using Solexa (Illumina) technology. Total 4,660,291 and 4,389,132 raw tags were generated in the control and BmCPV infected DGE libraries, respectively. A total of 4,498,263 and 4,258,240 clean tags corresponding to 49,365 and 60,834 distinct clean tags were filtered from the raw tags (Table 1). The distribution of the total and distinct tags over the different tag abundance categories showed highly similar patterns in each DGE library. More than 85% of the total clean tags had a copy number higher than 100 counts, whereas less than 7% of distinct tags had a copy number higher than 100 counts (Fig. 1). The sequencing saturation was analyzed to estimate whether the sequencing depth was sufficient for transcriptome coverage. The results showed that the number of detected genes was almost saturated when the total tag number reached 2 million or higher. Our sequencing depths reached approximately 3.5 million in each DGE library (Fig. S1), which satisfied the requirement for the experiment. So the two DGE libraries were reliable. 3.3. Analysis of tag mapping We mapped the tag sequences of the two DGE libraries to the reference database of silkworm, which contains 16,329 reference genes including 13,328 genes with the CATG site. In the control and BmCPV-infected DGE library, 31.92% and 29.61% of the clean tags corresponding to 20.3% and 19.71% of distinct tags were
10
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13
Table 1 Distributions of tags in two DGE libraries. Category
Parameter
Control (4008)
BmCPV-infected (4008)
Raw data
Total Distinct tags Total number Distinct tag number Total number Total% of clean tags Distinct tag number Distinct Tag% of clean Total number Total% of clean tags Distinct tag number Distinct Tag% of clean Number % of ref genes Number % of ref genes Total number Total% of clean tags Distinct tag number Distinct Tag% of clean Total number Total% of clean tags Distinct tag number Distinct Tag% of clean
4660291 122861 4498263 49365 1435964 31.92% 10023 20.3% 792073 17.61% 9387 19.02% 3956 24.23% 3728 22.83% 1547514 34.4% 17254 34.95% 840616 18.69% 14132 28.63%
4389132 144664 4258240 60834 1260841 29.61% 11991 19.71% 851711 20.00% 11372 18.69% 4621 28.3% 4359 26.69% 1739966 40.86% 22582 37.12% 755891 17.75% 16145 26.54%
Clean tags
All tags mapping to gene
Unambiguous tags Mapping to gene All tag-mapped genes Unambiguous Tag-mapped genes Mapping to genome
Unknown tags
tags
tags
tags
tags
Fig. 1. Distribution of total clean tags and distinct clean tags in each DGE library. C3: the control midgut of 4008 silkworm of the 3rd day in the fifth-instar; T3: the BmCPVinfected midgut of 4008 silkworm of the 3rd day in the fifth-instar. The numbers in square brackets indicate the range of copy numbers of each tag category. The data in parentheses indicate the number and percentage of corresponding tags among the total clean tags and distinct clean tags.
mapped to a gene in the reference database, 34.4% and 40.86% of the clean tags could be mapped to genome of silkworm, while 18.69% and 17.75% of the clean tags were unknown tags (Table 1). The number of unambiguous tags for each gene was normalized to TPM to accurately measure the gene expression level. The results showed that the mRNA transcribed from the major types of genes were represented in low copy number, and only a small proportion of genes were highly expressed (Fig. S2). 3.4. GO analysis and KEGG pathways A total of 752 differentially expressed genes were detected, of which 648 were upregulated and 103 were downregulated
(Fig. 2). To understand the functions of these differentially expressed genes, all of the 752 differentially expressed genes were mapped to terms in GO database and compared with the whole genome background. GO has three ontologies, namely, molecular function, cellular component, and biological process. 408 differentially expressed genes including 354 up-regulated genes and 54 down-regulated genes have a GO ID and can be categorized into a total of 41 functional groups using the WEGO software (Fig. 3). In each of the three main categories of the GO classification, ‘‘cellular process’’, ‘‘metabolic process’’, ‘‘cell and cell parts’’, ‘‘catalytic activity’’, ‘‘binding’’ terms are dominant. Interestingly, some genes related to the cell-proliferation, death, growth, locomotion, signaling, membrane-enclosed lumen, structural molecule, biological
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13
11
4986 significantly changed tags were found and mapped to 752 genes with 648 up-regulated and 103 down-regulated genes in DGE libraries. A total of 21 differentially expressed genes including 11 up-regulated and 10 down-regulated genes from DGE libraries were selected for qRT-PCR analysis to validate the DGE data. The results showed that 20 genes were demonstrated to have a concordant direction of change for both DGE and qRT-PCR except Bm_nscaf463_32 gene (Fig. 4).
4. Discussion
Fig. 2. Differentially expressed genes. The red part represents those genes upregulated in BmCPV-infected group compared to control group. The green part shows those genes down-regulated in BmCPV-infected group. The blue part shows those genes without expression difference between these two samples. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
process, nucleic acid binding, gene expression were all up-regulated in BmCPV-infected midgut of 4008 silkworm. KEGG (http://www.genome.jp/kegg) ontology assignments were used to classify the functional annotations of the identified genes to further understand the biological functions of the genes. Among the differentially expressed genes, 613 were mapped to 212 pathways in the KEGG database in 4008 DGE class. 19 terms was significantly enriched (Q value 60.05) (Table S2). 3.5. Identification and verification of differentially expressed genes The differentially expressed tags between the two samples were identified by FDR 60.001 and |log 2 ratio|P1. A total of
Bombyx mori cytplasmic polyhedrosis virus (BmCPV) is an important pathogen to the silkworm, always causing serious damage to sericulture production. However, the molecular mechanism of infection and defense response of the silkworm remains unclear. Silkworm strains p50 and 4008 are different in resistance to BmCPV infection according to previous studies. It showed that p50 is more resistant and 4008 is more susceptible (Xu et al., 2002; Wu et al., 2013). The previous microarray study by Wu et al. (2011) was focused on strain p50 but our present study was focused on the susceptible strain 4008. Two DGE libraries were made from midgut of the control and BmCPV-infected 4008 silkworm using Solexa (Illumina) technology with 4,660,291 and 4,389,132 raw tags, respectively. A total of 752 differentially expressed genes were detected, of which 648 were upregulated and 103 were downregulated (Fig. 2). The number of upregulated genes was markedly higher (over 6-fold) than the number of downregulated genes, which is different to the previous microarray study of B. mori larvae p50 strain infected with BmCPV with 135 upregulated and 123 downregulated genes in hybridization-based (microarrays) technology. In fact, we also examined the transcriptome profiling between the infected and mock-infected larvae of strain p50 with the same Digital Gene Expression approach and detected not only less number of differentially expressed genes but also less number of up-regulated genes than down-regulated ones, which is consistent to the result of the previous microarray study (data not provided). Therefore, the differential expression of genes in the two silkworm strains might be
Fig. 3. Histogram presentation of Gene Ontology classification. The results are summarized in three main categories: biological process, cellular component and molecular function. The y-axis on the left indicates the percentage of a specific category of genes. The y-axis on the right indicates the number of genes in a category.
12
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13 15 gene expression in qRT-PCR gene expression in DGE
Relative expression level (Fold change=log2 Ratio)
10 5
0
-5 -10
-15
Genes Fig. 4. Differential expression levels of candidate genes in BmCPV-infected and control midgut of 4008 silkworm larvae. The y-axis indicates the relative expression level of candidate gene mRNA transcripts (Fold change = log 2 Ratio). The x-axis indicates the candidate genes (from left to right) with the names of Bm_nscaf1690_036, Bm_nscaf2829_122, Bm_nscaf2964_026, Bm_nscaf2902_170, Bm_nscaf3031_071, Bm_nscaf2891_038, Bm_nscaf463_32, Bm_nscaf3062_05, Bm_nscaf2204_100, Bm_nscaf2529_090, Bm_nscaf2204_124, Bm_nscaf2829_149, Bm_nscaf2828_137, Bm_nscaf3075_32, Bm_nscaf98_23, Bm_nscaf2901_47, Bm_nscaf3099_033, Bm_nscaf3058_068, Bm_nscaf2575_003, Bm_nscaf2887_005 and Bm_nscaf2529_052. The genes with different expression regulation by qRT-PCR were indicated with diamonds ( ). Vertical bars represent the mean ± SE (n = 3), and bars in different time marked with star were significantly different from each other (p < 0.01).
involved in different resistance or defense response to BmCPV infection, which need to be clarified in our future research. Among the differentially expressed genes, 613 genes were mapped to 212 pathways in the KEGG database. It showed that the ribosome and RNA transport pathways were significantly enriched. Genes related to the terms ‘‘metabolic process,’’ ‘‘structural molecule activity,’’ and ‘‘cell or cell part’’ were dominant in the BmCPV-infected 4008 silkworm. Meanwhile, genes related to ribosome, gene expression, translation, and ribonucleoprotein complex biogenesis were all upregulated (Fig. 3). This result is mostly similar with the result of a DNA microarray employed by Wu et al. (2011). It suggested that BmCPV infection extensively affects known gene expression pathways. During the period of BmCPV infection, a large number of genes used for translation were overexpressed to meet the needs of the enormous proliferation of BmCPV. Among the 752 differentially expressed genes, 11 upregulated and 10 downregulated genes were verified via qRT-PCR (Table S1). Among the 20 concordant genes, only 12 genes, Bm_nscaf1690_036, Bm_nscaf2964_026, Bm_nscaf3031_071, Bm_nscaf3062_05, Bm_nscaf2204_100, Bm_nscaf2529_090, Bm_nscaf2829_149, Bm_nscaf98_23, Bm_nscaf2901_47, Bm_nscaf3058_068, Bm_nscaf2575_003, Bm_nscaf2887_005 were identified with FDR 60.001 and |log 2 ratio|P1 verified using qRT-PCR (Fig. 4). Lepidopteron species can resist baculovirus infection through selective apoptosis of infected cells in the midgut epithelium and through sloughing off infected cells (Clarke and Clem, 2002). Thus, the upregulated expression of the programmed cell death (PCD) protein 5 gene (Bm_nscaf3031_071) in the midgut cells may induce apoptosis in BmCPV-infected silkworm (Chen et al., 2013, 2001; Wang et al., 2007). Meanwhile, PCD is mediated via a nongenomic pathway that includes Ca(2+) as a second messenger and the activation of protein kinase C/caspase-3-like protease (Manaboon et al., 2009); Bm_nscaf1681_279, which encodes a Calreticulin (CRT) gene, is one of the major Ca(2+) binding chaperone proteins of the endoplasmic reticulum (ER) and plays a critical role in the activation of Ca(2+)-dependent apoptosis (Lim et al., 2008). In the cells over-expressing CRT, the intracellular calcium concentration was significantly increased and the activity of protein kinase c mRNA were reduced. It suggested that the up-regulated expression of CRT can disrupts intracellular calcium regulation leading to a calcium-dependent apoptosis in the BmCPV-infected cells. Furthermore, Cytosolic Ca(2+) concentration is crucially affected by the activity of the intracellular store (the sarcoplasmic reticulum, SR), which regulates Ca(2+) release (Berridge et al., 2003; McCarron et al., 2004). Two distinct classes of Ca(2+) release
channels were the inositol 1,4,5-trisphosphate (IP3) receptor and the ryanodine receptor. The important modulatory proteins that interact with these channels were FK506 binding proteins (FKBPs), including FKBP12 and its isoform-FKBP12.6 (Van Acker et al., 2004; Zissimopoulos et al., 2007). These regulatory proteins are reportedly important modulators of intracellular Ca(2+) release. In addition, Cyclic ADP-ribose (cADPR) binds to FKBP 12.6 on rat islet ryanodine receptor and that the binding of cADPR to FKBP12.6 frees the ryanodine receptor from FKBP12.6, causing it to release Ca(2+) (Noguchi et al., 1997). In addition, recent study have showed that the Ca(2+) entry channels in prothoracic gland cells of B. mori are probably modulated by protein kinases A and C. Application of protein kinases A and C inhibitors to the prothoracic glands cells of the silkworm, B. mori, resulted in slow and gradual increases in intracellular Ca(2+) (Dedos et al., 2008). FK506-binding protein (FKBP), which has also been identified as an inhibitor of protein kinase C (Nakazawa et al., 2005) may play a role in slow and gradual increases in intracellular Ca(2+). So the up-regulated expression of FK506-binding protein gene (Bm_nscaf1898_230) and protein kinase c inhibitor gene (Bm_nscaf3090_7) probably induced the increase of intracellular Ca(2+) in the BmCPV-infected cells and triggered the cell death. Early and pivotal events in apoptosis are now known to occur in the mitochondria and the endoplasmic reticulum(ER), and release of cytochrome c from the mitochondria and calcium from the endoplasmic recticulum into the cytosol is considered to be requisites for apoptosis in response to different stimuli (Wang and ElDeiry, 2004). FK506-binding protein, which has been identified as an inhibitor of protein kinase C is also an important modulatory protein that interacts with Ca(2+) release channels and may play a role in slow and gradual increases in intracellular Ca(2+). Calreticulin (CRT) is one of the major Ca2+ binding chaperone proteins of the ER and an unusual luminal ER protein. In the cells overexpressing CRT, the intracellular calcium concentration was significantly increased and the activity of PKC was reduced, leading to calcium-dependent apoptosis. The function of these two genes in the calcium-dependent apoptosis in the silkworm and furthermore the host response to BmCPV by the calcium-dependent apoptosis should be clarified in our future study. Genes related to activity of digestive enzyme were all down-regulated in the BmCPV-infected silkworms. Those genes include, for example, a chymotrypsin-like serine protease (Bm_nscaf2901_47), an acidic lipase (Bm_nscaf2829_149) and an lipase-1 (Bm_nscaf2529_052) (Fig. 4) etc. BmCPV is specific for epithelial cells in the silkworm midgut, which acts as the main site for virus
K. Gao et al. / Journal of Invertebrate Pathology 115 (2014) 8–13
multiplication (Miyajima and Kawase, 1968). As disease progresses, white wrinkles typically occur in the posterior part of the midgut, which may affect the digestive and the absorptive functions of the midgut. So the fatty acid binding protein gene (Bm_nscaf3058_068) and sugar transporter gene (Bm_nscaf2575_003) were also downregulated in the BmCPV-infected silkworm. In conclusion, we obtained 752 differentially expressed genes from the control and BmCPV-infected midgut of 4008 silkworm through DGE analysis. Analysis results of the KEGG pathway showed that the ribosome and RNA transport pathways were significantly enriched. The expression profiles of some important differentially expressed genes involved in signaling, gene expression, metabolic process, cell death, binding, and catalytic activity changed. However, most of these differentially expressed genes were up-regulated. These results suggest that the upregulated expression of CRT, FK506-binding protein, and protein kinase c inhibitor gene probably led to a calcium-dependent apoptosis in the BmCPV-infected cells. Therefore, the death or apoptosis program may be an antiviral response to defend the host from BmCPV proliferation. But further studies with other techniques are still needed, Such as over-expression, RNAi, western blotting, co-immunoprecipitation, and so on. The current results provide important clues for future research on the molecular mechanism of BmCPV infection. Acknowledgments This work was supported by the Jiangsu Postdoctoral Science Foundation of China (No. 1202026C) and the National Basic Research Program of China (Grant No. 2012CB114600). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jip.2013.10.016. References Benjamini, Y., Yekutieli, D., 2001. The control of the false discovery rate in multiple 348 testing under dependency. Ann. Stat. 29, 1165–1188. Berridge, M.J., Bootman, M.D., Roderick, H.L., 2003. Calcium signalling: dynamics, homeostasis and remodelling. Nat. Rev. Mol. Cell Biol. 4, 517–529. Chen, Y., Sun, R., Han, W., Zhang, Y., Song, Q., Di, C., Ma, D., 2001. Nuclear translocation of PDCD5 (TFAR19): an early signal for apoptosis? FEBS Lett. 509, 191–196. Chen, C.H., Jiang, Z., Yan, J.H., Yang, L., Wang, K., Chen, Y.Y., Han, J.Y., Zhang, J.H., Zhou, C.M., 2013. The Involvement of Programmed Cell Death 5 (PDCD5) in the regulation of apoptosis in cerebral ischemia/reperfusion injury. CNS Neurosci. Ther. 19, 566–576. Clarke, T.E., Clem, R.J., 2002. Lack of involvement of haemocytes in the establishment and spread of infection in Spodoptera frugiperda larvae
13
infected with the baculovirus Autographa californica M nucleopolyhedrovirus by intrahaemocoelic injection. J. Gen. Virol. 83, 1565–1572. Cloonan, N., Grimmond, S.M., 2008. Transcriptome content and dynamics at singlenucleotide resolution. Genome Biol. 9, 234. Dedos, S.G., Kaltofen, S., Birkenbeil, H., 2008. Protein kinase A and C are ‘‘Gatekeepers’’ of capacitative Ca2+ entry in the prothoracic gland cells of the silkworm, Bombyx mori. J. Insect Physiol. 54, 878–882. Hill, C.L., Booth, T.F., Stuart, D.I., Mertens, P.P., 1999. Lipofectin increases the specific activity of cypovirus particles for cultured insect cells. J. Virol. Methods 78, 177– 189. Livad, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25, 402–408. Lim, S., Chang, W., Lee, B.K., Song, H., Hong, J.H., Lee, S., Song, B.W., Kim, H.J., Cha, M.J., Jang, Y., Chung, N., Choi, S.Y., Hwang, K.C., 2008. Enhanced calreticulin expression promotes calcium-dependent apoptosis in postnatal cardiomyocytes. Mol. Cells 25, 390–396. Magnoler, A., 1974. Effects of a cytoplasmic polyhedrosis on larval and postlarval stages of the gypsy moth, Porthetria dispar. J. Invertebr. Pathol. 23, 263–274. Manaboon, M., Iga, M., Iwami, M., Sakurai, S., 2009. Intracellular mobilization of Ca2+ by the insect steroid hormone 20-hydroxyecdysone during programmed cell death in silkworm anterior silk glands. J. Insect Physiol. 55, 122–128. McCarron, J.G., Bradley, K.N., MacMillan, D., Chalmers, S., Muir, T.C., 2004. The sarcoplasmic reticulum, Ca2+ trapping, and wave mechanisms in smooth muscle. News Physiol. Sci. 19, 138–147. Miyajima, S., Kawase, S., 1968. Changes in virus-infectivity titer in the hemolymph and midgut during the course of a cytoplasmic polyhedrosis in the silkworm. J. Invertebr. Pathol. 12, 329–334. Morozova, O., Marra, M.A., 2008. Applications of next-generation sequencing technologies in functional genomics. Genomics 92, 255–264. Nakazawa, T., Takasawa, S., Noguchi, N., Nata, K., Tohgo, A., Mori, M., Nakagawara, K., Akiyama, T., Ikeda, T., Yamauchi, A., Takahashi, I., Yoshikawa, T., Okamoto, H., 2005. Genomic organization, chromosomal localization, and promoter of human gene for FK506-binding protein 12.6. Gene 360, 55–64. Noguchi, N., Takasawa, S., Nata, K., Tohgo, A., Kato, I., Ikehata, F., Yonekura, H., Okamoto, H., 1997. Cyclic ADP-ribose binds to FK506-binding protein 12.6 to release Ca2+ from islet microsomes. J. Biol. Chem. 272, 3133–3136. Van Acker, K., Bultynck, G., Rossi, D., Sorrentino, V., Boens, N., Missiaen, L., De Smedt, H., Parys, J.B., Callewaert, G., 2004. The 12 kDa FK506-binding protein, FKBP12, modulates the Ca(2+)-flux properties of the type-3 ryanodine receptor. J. Cell Sci. 117, 1129–1137. Wang, S., El-Deiry, W.S., 2004. Cytochrome c: a crosslink between the mitochondria and the endoplasmic reticulum in calcium-dependent apoptosis. Cancer Biol. Ther. 3, 44–46. Wang, N., Lu, H.S., Guan, Z.P., Sun, T.Z., Chen, Y.Y., Ruan, G.R., Chen, Z.K., Jiang, J., Bai, C.J., 2007. Involvement of PDCD5 in the regulation of apoptosis in fibroblast-like synoviocytes of rheumatoid arthritis. Apoptosis 12, 1433–1441. Wang, Z., Gerstein, M., Snyder, M., 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57–63. Wu, P., Qian, H.Y., Xu, A.Y., Chen, T., Hou, Q.R., Guo, X.J., Li, L., 2013. Analysis on resistance of various Bombyx mori Varieties to cytoplasmic polyhedrosis virus infection by quantitative PCR. Sci. Sericul. 39, 0289–0294. Wu, P., Wang, X., Qin, G.X., Liu, T., Jiang, Y.F., Li, M.W., Guo, X.J., 2011. Microarray analysis of the gene expression profile in the midgut of silkworm infected with cytoplasmic polyhedrosis virus. Mol. Biol. Rep. 38, 333–341. Xu, A.Y., Li, M.W., Lin, C.Q., Zhang, Y.H., Hou, C.X., 2002. The comparision of the resistibility to Bombyx mori cytoplasmic polyhedrosis virus among silkworm germplasm resourse. Sci. Sericul. 28, 157–159. Ye, J., Fang, L., Zheng, H., Zhang, Y., Chen, J., Zhang, Z., Wang, J., Li, S., Li, R., Bolund, L., 2006. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 34, W293–W297. Zissimopoulos, S., Docrat, N., Lai, F.A., 2007. Redox sensitivity of the ryanodine receptor interaction with FK506-binding protein. J. Biol. Chem. 282, 6976–6983.
