Baculovirus Induced Transcripts in Hemocytes from the Larvae of Heliothis virescens

Using RNA-seq digital difference expression profiling methods, we have assessed the gene expression profiles of hemocytes harvested from Heliothis virescens that were challenged with Helicoverpa zea single nucleopolyhedrovirus (HzSNPV). A reference transcriptome of hemocyte-expressed transcripts was assembled from 202 million 42-base tags by combining the sequence data of all samples, and the assembled sequences were then subject to BLASTx analysis to determine gene identities. We used the fully sequenced HzSNPV reference genome to align 477,264 Illumina sequence tags from infected hemocytes in order to document expression of HzSNPV genes at early points during infection. A comparison of expression profiles of control insects to those lethally infected with HzSNPV revealed differential expression of key cellular stress response genes and genes involved in lipid metabolism. Transcriptional regulation of specific insect hormones in baculovirus-infected insects was also altered. A number of transcripts bearing homology to retroviral elements that were detected add to a growing body of evidence for extensive invasion of errantiviruses into the insect genome. Using this method, we completed the first and most comprehensive gene expression survey of both baculoviral infection and host immune defense in lepidopteran larvae.


Introduction
Baculoviruses are insect-specific pathogens that have been utilized as both tools for recombinant protein expression and biological agents targeting agriculturally significant pests, most importantly within the order Lepidoptera [1][2][3][4]. Helicoverpa zea single nucleopolyhedrovirus (HzSNPV) infects several species of noctuid moths with outcomes ranging from a low rate of infection to 100% mortality, depending on the host species [5][6][7]. This range of outcomes points to the importance of the host response to viral infection, about which little is understood at the molecular level, within the infected animal. Methodologies that employ DNA microarray and gene chip have been utilized to probe pathogen-host interactions at the transcriptional level to shed light on factors that determine the outcome of infection, but these methods rely on existing consensus DNA sequence information for a given species [8][9][10]. While such studies are possible in genetically well-defined experimental models, this limitation has restricted genetic and transcriptional analysis of less intensively studied hosts and their interactions with relevant pathogens [11][12][13]. Indeed, the bulk of insect host response to viral infection has been related to the dipteran response to infection with RNA viruses, chiefly Drosophila melanogaster [14][15][16][17] and mosquitoes [18][19][20].
RNA-seq studies of host gene expression allow quantitation of gene expression across an entire genome in response to infection [21,22]. Examination of host gene expression in response to baculovirus infection has been the subject of relatively few studies [23]. Based on our prior observations, and those of others, we hypothesized that baculovirus infection of a permissive species would result in substantial changes in the transcriptional profile of the host. To address this hypothesis in a biologically relevant in vivo context, both baculovirus and host gene expression were documented simultaneously in the same tissues of infected Heliothis virescens larvae by an RNA-seq approach using an Illumina GAII sequencing instrument. Illumina is a sequencing method that generates up to a terabase of sequence information as short (approximately 40 bases) sequence reads from input samples of nucleic acids [23]. By barcoding the nucleic acids of each sample, multiple treatments, time points, controls and different species were crowded onto a single lane, allowing a fine grained understanding of the regulation of thousands of genes simultaneously. Illumina RNA-seq has been used to examine differential gene expression during whitefly infection by tomato yellow leaf curl china virus, and the corn planthopper subsequent to infection with Maize mosaic rhabdovirus [24,25].
After assessing sample quality and screening out contaminant rRNA and tRNA, we constructed a de novo transcriptome of approximately 100,000 Heliothis virescens hemocyte transcripts and subsequently aligned sequence tags from control and HzSNPV-infected hemocyte treatments to this reference transcriptome. Using a digital expression profiling approach comparing control and infected samples, we observed overall trends in important cellular processes closely related to the natural history of infection, and expanded on previous studies that focused more narrowly on the differential expression of specific host genes consequent to infection [26]. Taken together, this study surveys the comparative transcriptional state of infected larvae to reveal pathogen-specific alterations in gene expression and simultaneously represents the largest deposition of H. virescens mRNA sequence information in a publicly-accessible database to date.

Transcriptome Generation and Assembly
The treatment and enumeration of the Illumina sequence output is summarized in Table 1. A total of 202 M 32-42 nt base tags from each of the three treatments and control lanes (Table 1, Total Seqs) were pooled, cleaned, and screened (Table 1, After Decontamination) for tRNA and rRNA sequences and assembled as described in the Experimental Section. A reference de novo transcriptome assembly of 103,879 sequences was thus generated from this output of hemocytes from HzSNPV infected, bacterial-and fungal-stimulated larval hemocytes. Approximately 31 million of the control reads aligned to this reference assembly, and of these almost 20 million aligned uniquely. While mRNA derived from hemocytes of bacterial and fungal septic puncture treated insects were utilized in constructing the de novo transcriptome assembly, the subsequent gene regulation data for those treatments are outside the scope of this manuscript and will be reported elsewhere.

Illumina Reads Matching the HzSNPV Genome
Of the 202 million reads used to assemble the Illumina library, 477,264 aligned with the HzSNPV genome (approximately 0.26%) [27]. The majority of the reads aligned uniquely within one of the 139 HzSNPV open reading frames (ORF) ( Figures 1A,B). The ORF with the highest number of matching reads was ORF31 38K/pp31 (34,047 matches) followed by ORF11 odv-ec27 (22,096 matches) which correspond to the Autographa californica multiple nucleopolyhedrovirus (AcMNPV) ORFs 36 and 144, respectively ( Figure 1A) [28]. Many of the ORFs with the lowest number of matching reads were ORFS that are to date specific to HzSNPV and/or Helicoverpa armigera SNPV ( Figure 1B) [27]. Only 4064 Illumina reads aligned to more than one ORF (Table 2), and these reads were spread between 14 ORFs. The reads that aligned to more than one ORF did so because of overlapping areas between two ORFs. The difference in number of hits refers to the difference between the number of reads that are unique to the ORF and those that match to the area of two overlapped ORFs.
Four ORFS have previously been found that are unique to HzSNPV (ORF26, ORF42, ORF62, and ORf79) [27]. ORF79 had two Illumina read matches which was the lowest number of reads for any ORF and suggests that it is not a functional ORF. ORF26 had 29 hits but seven of these potentially overlapped ORF27. ORF62 had 74 hits, which was greater than the number of matches for ORFs 22, 107 and 15. ORF42 had 354 matching reads, which was more than lef-10 (298 hits) and encodes for a similar size protein, adding to the likelihood that it is a functional ORF.

HSP70/Small HSPs
In hemocytes isolated from insects infected with baculovirus, transcription of heat-shock protein 70 (hsp70) genes was up-regulated 70-and 176-fold over control levels, and two transcripts comprised of 1763 and 1452 nt were assembled that bore homology to hsp70 of Spodoptera exigua and Helicoverpa zea, representing potentially different isoforms of this gene (Table 3). Four small HSP transcripts encoding HSPs of 19 to 21 kiloDaltons (kD) were identified as upregulated between 16-and approximately 21-fold (Table 3).

Cellular Adhesion/Immunity Genes
Baculovirus infection suppressed transcription of cell adhesion molecules dystroglycan, and fasciclin-1 to 0.3 and 0.02 (three-and 50-fold) relative to control levels, respectively, whereas a gene identified as EGF-containing fibulin-like extracellular matrix protein 2 was induced seven-fold over control ( Table 3). Transcription of the immunity-related signal transducing adapter molecule 2 was suppressed approximately 4.5-fold, while cecropin D was upregulated by 12-fold in response to viral infection. Larval cuticle protein 2 was also up-regulated following baculovirus infection.

Endocrine Transcripts
A gene identified as probable nuclear hormone receptor HR3 was upregulated nearly eight-fold in response to baculovirus infection (Table 4). Hormone receptor 3C, juvenile hormone epoxide hydrolase, and nuclear hormone receptor FTZ-H1F were upregulated between approximately three-to four-fold. Transcription of basic juvenile hormone-suppressible protein 1 and 2 were suppressed in the presence of virus by 40-and 200-fold, respectively; juvenile hormone diol kinase was also down-regulated following baculovirus infection. Table 3. Differential regulation of stress response, adhesion, and immunity-related genes.

Metabolic Regulation
mRNA sequences bearing significant homology to stearoyl-CoA desaturase, fatty acid reductase, and fatty acid synthase genes were upregulated by baculovirus infection (Table 4). A total of sixteen genes related to lipid metabolism were induced by viral infection with magnitudes of induction in response to baculovirus that ranged from 3-to 237-fold. Four genes encoding enzymes involved in this pathway were down-regulated in response to baculovirus infection; among these was fatty acid binding protein 3, which functions to bind to unsaturated fatty acids. Genes identified as p260/p270, thought to be fatty acid synthase homologs, were upregulated in response to HzSNPV infection: p260 by more than 400-fold, and p270 by 51-fold. Related to this metabolic pathway, the transcription of genes encoding storage proteins arylphorin, p82 (riboflavin-binding hexamer), and hexamerine was suppressed between 20-and 100-fold.

Errantivirus
More than a hundred genes that were differentially regulated in response to baculovirus infection were identified as homologs of endogenous retroviral elements or retrotransposon related genes. In insects, these are termed errantiviruses [29]. Because of the large number of these transcripts, an arbitrary reporting requirement of at least 45 sequence detections of a given gene in either the control or virus-infected samples was established. With one exception, the remaining genes were upregulated between 3.3-and 11.25-fold compared to those detected in the control hemocytes (Table 5).

Insects, Infection, and RNA Isolation
Heliothis virescens eggs were received from the North Carolina State University, Dept. of Entomology Insectary (Raleigh, NC, USA). Larvae were reared individually on an artificial wheat germ based diet (Bio-Serv, Frenchtown, NJ, USA) under standard conditions of 14 h:10 h (L:D) photoperiod, 55% RH, 28 °C [30]. Synchronized early-4th instar larvae (4 hours post molt) were either mock-infected with a control inoculum absent virus or infected with an LC99 dose of HzSNPV per os (5 × 10 8 polyhedra/mL) [30]. To activate the antibacterial immune response early 5th instar larvae were punctured with a tungsten needle dipped into a 1 μg/mL suspension of lipopolysaccharide and peptidoglycan (Sigma Chem. Co., St. Louis, MO, USA) in PBS [31]. The antifungal response was activated by puncture with a tungsten needle dipped into a 1 μg/mL suspension of β-glucan, curdlan and laminarin (Sigma Chem. Co., St. Louis, MO, USA) in PBS. Hemocytes were collected at 12 hpi in order to capture early host-responses to HzSNPV infection. Hemocytes were collected from 30 insects from each treatment according to previously published methods [32]. Larvae were bled through a punctured anterior proleg into ice cold PBS containing a crystal of phenylthiourea to prevent melanization. Hemocytes were pelleted by centrifugation at 5000 × g for 4 min, plasma supernatant was removed, and stored at −85 °C for later use. Total RNA was extracted from hemocytes using RNeasy™ kits (Qiagen, Valencia, CA, USA) [33].

Sequence Generation, Assembly and Annotation
Hemocyte RNA pools were submitted to University of Missouri, Bond Life Sciences Center DNA Core for Illumina GAII sequencing. RNA quality of the pooled time points was examined using the Experion Automated Electrophoresis system (Bio-Rad, Hercules, CA, USA). Libraries were constructed according to the standard Illumina RNA-seq protocol (Part# 1004898 Rev. A, rev Sept 08) for the pooled PCR products except for the fragmentation step [34]. First strand and second strand cDNA synthesis was carried out with Invitrogen superscript II reverse transcriptase and random hexamer primers according to the manufacturers protocol. Illumina paired-end adapters were ligated and a gel size-selected ~300 bp fragments depleted of adapter dimer and large chimeric fragments were subjected to 15 rounds of selective PCR enrichment. The four treatments were run on four separate lanes on an Illumina GAIIx sequencing instrument. To assemble the de novo hemocyte transcriptome the sequence data from all lanes was combined yielding over 202 million 42-base long Illumina single end reads which were first cleaned of low quality bases at the 3' end, and then cleaned for low quality overall using the Fastx Toolkit [35]. The remaining reads were screened for homology to contaminants (predominantly rRNA) to leave 192 million reads of lengths 32 to 42. VELVET was used to assemble the reads into preliminary contigs using a range of word size (K values) from 21 to 41 [36]. OASIS was used to assemble the VELVET output into transcript isoforms of at least 100 bases in length [37]. Those reads that were not assembled under the first round were used as input to a second round of assembly by VELVET and OASIS. The results for each value of the word size were combined and redundant transcripts were removed using Vmatch resulting in 103,879 unique potential transcripts [38]. These sequences were annotated by homology using BLASTx with proteins in the UniProt uniref100 database [39]. About two dozen sequences were identified as chimeric and removed. Functional categories of contig and singleton sequences resulting from this assembly were annotated using a local installation of BLAST2GO [40]. All sequence data discussed in this manuscript have been deposited in NCBI GenBank [41].

Alignment of Sequence Tags with HzSNPV Reference Genome
The cleaned reads from the various samples were aligned to the HzSNPV genome using SOAP [42]. ORFs were defined from reported coordinates at GenBank [27].

Identification of Differentially Regulated Genes
Host gene expression was analyzed in a semi-quantitative fashion by comparison of the control and HzSNPV Illumina lanes, normalizing the number of times sequence tags aligning with an mRNA was detected in the cells of HzSNPV infected larvae to those of controls. We analyzed genes that were shown to be upregulated or down regulated by a factor of three and that had an e-value less than 0.0001, and those findings are presented in Tables 3 through 5.

Discussion and Conclusions
Previous work has shown that during HzSNPV infection of H. virescens larvae, infection of midgut columnar epithelial cells could be detected by 4 hours post infection but remained constant by 24 hours post infection [43,44]. In larvae infected with HzSNPV-hsp70/LacZ using a dose killing 84% of a test group, 40% of the larvae were LacZ positive in either the midgut or midgut and tracheal epithelia [43]. In this study, we applied an RNA-seq approach to examine transcription of hemocytes extracted from H. virescens that were experimentally challenged with HzSNPV to capture the early transcriptional profile of the baculovirus genome within the first 12 hours of infection ( Figures 1A,B). At 12 hours post ingestion, we expected to capture a time point at which the primary infection of the midgut had begun to progress into the circulating hemocytes, in attempt to identify the host genes that are potentially involved in an immune response to the virus. Surprisingly, we detected transcription of all known HzSNPV ORFs at this time point at varying levels of expression. While this analysis provides some detail of the level to which these genes are transcribed early in infected hemocytes, more infection time points would be necessary to precisely follow gene expression levels and allow a full understanding of the temporal interaction of host and viral genes during baculovirus infection.
A comparison of mRNA extracted from HzSNPV-infected larvae and those of mock-infected larvae revealed a number of differentially regulated genes of physiological and immunological relevance. Induction of the heat shock protein 70 gene (hsp70) in response to productive baculovirus infection was previously demonstrated by our lab using 2D-gel electrophoresis analysis and quantitative PCR studies of infected lepidopteran insects and lepidopteran cell lines [26,45]. Our finding of hsp70 induction in response to baculovirus infection also correlates with recent studies indicating that HSP70 plays an important role in baculovirus expression vector systems [46]. This study shows that virus infection results in the induction of hsp70 of up to a 177-fold. The heat shock protein family of proteins are well-conserved across the genera and contribute to proper protein folding, organization, and function, especially during stress induced by thermal, chemical, or pathogen-induced cellular perturbations. The hypothesis that baculovirus exploits HSP70 to enhance the availability of properly folded viral proteins during its replication cycle is currently under further study.
In mammals, adhesion modules play a number of roles in response to pathogens, at times serving as a host-defense mechanism, and in other cases they are exploited by pathogens to facilitate dissemination within the host or dispersal without [47][48][49]. Substantial changes were detected in the expression of cell adhesion molecules fasciclin-1 and dystroglycan, as well as a probable matrix anchor protein, following baculovirus infection. Lepidopteran homologs of fasciclin-1 and dystroglycan have not been well characterized, and what role they play during baculovirus infection remains speculative. Homologs of fasciclin-2 and fasciclin-3, which have roles in neural and midgut development [50][51][52][53] were detected in our study, but the level of their respective transcripts was not altered as compared to controls in any of the hemocytes from the challenged insects (data not shown).
Cecropin D, an important molecule for bacterial immunity, was also induced in our study by baculovirus infection [54,55]. A suppression of signal transducing adapter molecule 2 expression of approximately five-fold was noted, although transcription of the signal transducing activator of transcription (stat) gene, which plays a critical role in animal immune signaling, was unperturbed (data not shown) in response to any of the treatments.
Fatty acid metabolism and hormone regulation was also altered in response to baculovirus infection. The fatty acid synthase homologs p260/p270 play a role in limb development in Bombyx mori [56], and were highly upregulated in response to baculovirus, but not in response to either fungus or bacteria injection. This trend of up-regulation of fatty acid metabolism-related transcripts was generally reflected among the sampled genes, with a number of lipid reductases and desaturases upregulated more than 30-fold following baculovirus infection. Also uniquely reported here, but perhaps not surprisingly, was the magnitude and uniformity of the observed down-regulation of storage proteins hexamerin/arylphorin in response to baculovirus challenge. Commensurately, transcription of a gene identified as p82 riboflavin-binding hexamer monomer was suppressed 20-fold in response to baculovirus infection. A role for the extensive modulation of fatty acid metabolism in the replication of baculovirus has yet to be determined, but this apparent shunting of energetic resources from a state of growth to one of infection response has previously been reported in polydnavirus-infected H. virescens as well as Serratia marcescens-infected honey bees [57,58].
Alterations in transcript levels of genes related to hormonal activity demonstrated the close association between baculovirus infection and the developmental hormones in larvae. A transcript exhibiting homology to the nuclear hormone receptor HR3, inducible by ecdysteroid, was upregulated following baculovirus infection. Baculoviruses are known to decrease the ecdysteroid levels in larval hemolymph by glycosylating ecdysteroid through the ecdysteroid UDP-glycosyltransferase (EGT) gene thereby causing larvae to delay pupation [59]. Juvenile hormone (JH) epoxide hydrolase and basic JH suppressible homologs were inversely regulated: transcription of the former increased between three-and four-fold, while transcription of the latter was suppressed up to 200-fold. Decreasing the JH levels in larvae signals the cessation of feeding and the beginning of the molt process. Alterations or fine tuning of the enzyme levels involved in the JH removal would also increase/maintain the JH titer and again delay pupation allowing the baculovirus to continue to propagate. The overexpression of JH esterase in recombinant baculoviruses has been used to successfully improve their efficacy [60,61].
Errantiviruses have been detected in Drosophila, and in lepidopteran species B. mori and Lymantria dispar. Retroviral sequences were upregulated by baculovirus infection in cultured Sf9 insect cell lines utilized in BEV systems, derived many years ago from Spodoptera frugiperda, and also the Hi-5 cell line isolated originally from the cabbage looper Trichoplusia ni [62]. Numerous transcripts containing errantivirus and retrovirus-like sequences were detected in our study. It has been observed that insect retroviruses most likely derived portions of their genome from baculoviruses via recombination during their evolution [63]. Menzel and Rohrman speculate a plausible scenario in which baculovirus infection may trigger transposition events, either by a direct mechanism intrinsic to baculovirus(es) or by mitigating the suppressive effect of host silencing factors [62]. Indeed, errantiviruses have been known to "piggy back" by integrating into the baculovirus genome during infection [64], and our report provides further evidence for the activation of errantivirus during productive baculovirus replication in an infected insect.
Overall, this study demonstrates the utility of next-generation sequencing in characterizing host response to various treatments at the transcriptional level. In addition to capturing genetic sequences from a species lacking extensive sequence information, comparison between control and HzSNPV-infected insects has demonstrated key physiological differences in the transcriptional profile of the infected host. Trends were identified that point towards a decrease in host energy storage, alterations in transcription of hormonal genes affecting molting and pupation, altered transcription of genes related to cell adhesion and motility, and a generalized stress response that resulted in the virus-specific transcriptional induction of discrete heat shock response constituents. Analysis of baculovirus-specific transcripts has revealed the coordinated expression of key viral genes following infection, as well as a surprising number of late genes that appeared to be transcribed earlier than expected. Further, assembly of the contigs isolated from the hemocytes of baculovirus-infected insects netted a number of differentially regulated retroviral-transcripts that were deduced to be errantiviruses, of which very little is currently known. Further mining of these data and subsequent experimentation is expected to add to the depth and breadth of our understanding of the complex interactions that govern the outcome of baculovirus infection in their lepidopteran hosts.