Novel insights into the insect trancriptome response to a natural DNA virus

Larval culture, RNA extraction, library preparation and sequencing

The host species used in our experiment was Plodia interpunctella, the Indian meal moth. All individuals were taken from a large outbred stock that has been maintained at the University of Sheffield for approximately 8 years and reared on a cereal-based diet under standard laboratory conditions (27°C; 16 L:8D cycle). The pathogen used was Plodia interpunctella Granulosis Virus (PiGV), a natural DNA virus of P. interpunctella. The virus is naturally transmitted through the ingestion of virus occlusion bodies either from the environment or through necropsy of infected individuals. Baculoviruses (BVs) are generally host-specific, obligate killing DNA viruses [26] with two phenotypes; the occluded virus, which consists of virions encased within a protein coat, and the budded virus. The initial site of virus exposure in the insect is the midgut where the alkaline conditions destroy the occlusion body resulting in the release of infectious virions. These virions may cross the peritrophic membrane into the epithelial cells and systemic infection is established when progeny BVs pass across the basal lamina into the haemocoel and infect secondary tissue such as the tracheal system and the fat body [27]. Our aim was to capture the change in mRNA levels during early infection, as immune responses can be swift and transient [28]. To this end, we chose to sample at 24 hours post infection, at which point enveloped virus nucleocapsids were found trapped in the connective tissue surrounding the midgut [29].

Transcriptome assemblies were constructed from 18 P. interpunctella samples sequenced with RNA-Seq (Illumina) 100 base, pair-end reads on three HiSeq lanes (Illumina). The 18 samples came from two independent experiments, one of which (N = 12) is described in a separate paper (McTaggart et al., in prep). Each sample consisted of 20 pooled individuals. Here, our aim was to characterize the transcriptomes of larvae, 24 hours after they had been exposed to viral dose of LD₅₀. We set up six replicates in each of two treatments (viral exposed and control). In each replicate, 30 F1 generation adult moths were kept together with excess food. Prior to their viral exposure, 50 fourth instar larvae were removed from the pot, starved for 2 hours and then inoculated with the virus using a standard droplet feeding method [27]. The virus solution was prepared by centrifugation of infected P. interpunctella larvae, and diluted in blue food colouring and sucrose. The required dose was calculated from a dose response assay performed prior to the experiment.

We inoculated larvae with either a solution of purified PiGV to LD₅₀ or with a control solution consisting of only sucrose and blue food colouring. Both treatments were orally administered with a standard oral droplet feeding protocol, and considered successful if the blue solution was visible in at least half of the gut length. Inoculated larvae were moved to individual wells of a 25-cell Petri dish containing ample food. Twenty-four hours later twenty larvae from each replicate were pooled and crushed in 1 ml of Trizol (Life Technologies) and stored at −20°C until RNA extraction. The remaining larvae (N ~ 30) were monitored for infection for 30 days. Due to financial constraints, only three of the six replicates from each treatment (i.e. six samples in total) were sequenced using RNA-Seq.

Immediately before RNA extraction, an additional 500ul of TRIzol (Ambion, Life Technologies) was added to each sample, and incubated for 5 minutes at room temperature. Two hundred ul of chloroform was added to each sample, shaken vigorously for 15–20 seconds and then centrifuged at 11600 rcf for 15 minutes at 4°C. The upper, aqueous phase was isolated and nucleic acids precipitated by adding 0.5 volumes of Absolute Ethanol, and inverting the tubes several times. This solution was used as the starting material for the RNAeasy (Qiagen) protocol, which was followed according to the manufacturer’s instructions. The integrity of the resultant total RNA was confirmed on Bioanalyzer (Agilent RNA-nano reagents). RNA and DNA concentrations were determined with a Qubit fluorometer (Invitrogen QuantRNA), while the 260:280 ratio was assessed on a Nanodrop (ThermoScientific). For each sample, we subjected 5 ug of total RNA to one round of poly-A selection on oligo(dT) Serabeads. The resultant messenger RNA was fragmented to an average size of 100 bp using divalent cations at 95°C for 5 min prepared following the manufacturer’s recommended protocol (Illumina mRNAseq kits Cat no. RS-100-0801). First strand cDNA synthesis was carried out using Superscript III reverse transcriptase (Invitrogen) and 3 ug random hexamer primers (Illumina) per sample as per the manufactures’ instructions. Second strand cDNA synthesis and RNAseq samples were prepared according to the manufacturer’s recommended protocol (Illumina). The fragment size and concentration of resultant libraries were assessed on a Qubit fluorometer (Invitrogen QuantRNA) and on a Bioanalyser High Sensitivity Chip (Invitrogen QuantRNA).

Transcriptome assemblies

Adapter sequences were trimmed from the raw reads using the program Scythe (https://github.com/vsbuffalo/scythe). The program Sickle (https://github.com/najoshi/sickle) was used to remove low quality bases, reads with N’s, and sequences that were less than 50 bases long. SOAPdenovo-Trans [7] was used to build two assemblies with different k-mer sizes (kmer = 25, kmer = 30). Otherwise, default parameter settings were used. The assembler Trinity [8] was also run on the data twice, once with an edge threshold parameter of 0.05 and the other of 0.16. Default settings were used for all of the other parameters. Reads were mapped back to each assembly using GSNAP (http://research-pub.gene.com/gmap/). Reads thought to be due to PCR duplication were flagged with Picard (http://broadinstitute.github.io/picard/) and were not considered in any analyses. All contigs that had three or fewer reads mapping to it across all treatments were considered to be misassembled and removed.

In order to determine if any of the assembled contigs were of viral origin, the complete genome of Plutella xylostella granulovirus (NC_002593.1)(PxGV) was queried against the P. interpunctella transcriptome. Any contig that had a BLAST hit with an e-value of less than 1e-30 was considered to be viral. A dot-plot of these transcripts was made against the PxGV genome to estimate how much of the PiGV transcriptome had been recovered.

Assembly assessment

For the purposes of comparison, we assessed three basic parameters of each of the four assemblies: contig length, the number of contigs and the number of bases in the contigs. We also quantified the proportion of the sequence reads that mapped to each assembly using GSNAP. To test the completeness of the assemblies, we tested for the presence of highly conserved genes from other species in two ways. First, we built a BLAST database of each transcriptome that was queried with a previously published set of ultraconserved single-copy orthologs (UCOs) using tBLASTn. We counted how many times each gene was observed (with the expectation that each gene should be present only once) in each transcriptome. Similarly, we used the program CEGMA [30], which uses a different set of highly conserved single copy orthologs to query the data sets. The preferred assembly was then run through the Evidence Gene pipeline (http://arthropods.eugenes.org/genes2/about/EvidentialGene_trassembly_pipe.html). This pipeline reduces redundancy in de novo assemblies by translating each contig in all six reading frames, selecting the longest open reading frame. Pairwise comparisons between all contigs eliminate all that have the same protein sequences.

Functional annotation

We used the afterParty (https://github.com/mojones/AfterParty2) interface to annotate all of our contigs using BLAST to the Uniref90 database [31]. Canonical immune system genes from the RNAi, Toll, IMD and JNK pathways were manually annotated within, and are accessible at, the afterParty P. interpunctella database (http://afterparty.bio.ed.ac.uk/study/show/2194070). Contigs with a BLAST hit to the query having an e-value hit greater than 1e-10, and covering at least 80% of the full-length transcripts (inferred based on homology) were considered to be valid candidate transcripts. All contigs were loaded into Blast2Go, annotated with the non-redundant database of NCBI, and GO terms were assigned when possible. Combined graphs of the entire transcriptome were constructed for Biological Process and Molecular Function (level 2). These results were compared with the two other lepidopteran species for which equivalent data was available, Maruca vitrata [11] and Manduca sexta, with a chi-square test [12].

Differential expression of immune system genes

The number of times each contig in the chosen transcriptome was observed in the sequence reads was calculated using HTSeq (http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html). Differential expression between P. interpunctella that were exposed to PiGV compared to controls was calculated using DESeq (version 1.9.4) [32]. All genes (contigs) with an FDR correction of less than 0.10 were considered to be significant.

q-RT-PCR

We analysed the expression of comp6230_c0_seq1 (cuticle protein) and comp2004_c0_seq1 (GNBP) using comparative C_T (ΔΔC_T) qPCR with actin as the internal control gene. Using the RNA samples prepared for RNASeq, we synthesised complementary DNA (cDNA) by first mixing 500 ng total RNA with 500 ng random hexadeoxynucleotides (Promega) and heating to 75°C for 10 minute. The samples were chilled on ice after which 200 u MMLV reverse transcriptase (Promega), 5 μl 5x MMLV reverse transcriptase buffer, 1.25 μl dNTPs (Promega; final concentration 0.5 mM) and 20 u RNAsin RNase inhibitor (Promega) was added and the volume adjusted to 25 μl with nuclease free water. Samples were incubated at 37°C for 60 minutes, 70°C for 15 minutes and then stored at −20°C.