Transcriptome sequencing of and microarray development for a Helicoverpa zea cell line to investigate in vitro insect cell-baculovirus interactions

Abstract

The Heliothine insect complex contains some of the most destructive pests of agricultural crops worldwide, including the closely related Helicoverpa zea and H. armigera. Biological control using baculoviruses is practiced at a moderate level worldwide. In order to enable more wide spread use, a better understanding of cell-virus interactions is required. While many baculoviruses have been sequenced, none of the Heliothine insect genomes have been available. In this study, we sequenced, assembled and functionally annotated 29,586 transcripts from cultured H. zea cells using Illumina 100 bps and paired-end transcriptome sequencing (RNA-seq). The transcript sequences had high assembly coverage (64.5 times). 23,401 sequences had putative protein functions, and over 13,000 sequences had high similarities to available sequences in other insect species. The sequence database was estimated to cover at least 85% of all H. zea genes. The sequences were used to construct a microarray, which was evaluated on the infection of H. zea cells with H. Armigera single-capsid nucleopolyhedrovirus (HearNPV). The analysis revealed that up-regulation of apoptosis genes is the main cellular response in the early infection phase (18 hours post infection), while genes linked to four major immunological signalling pathways (Toll, IMD, Jak-STAT and JNK) were down-regulated. Only small changes (generally downwards) were observed for central carbon metabolism. The transcriptome and microarray platform developed in this study represent a greatly expanded resource base for H. zea insect-HearNPV interaction studies, in which key cellular pathways such as those for metabolism, immune response, transcription and replication have been identified. This resource will be used to develop better cell culture-based virus production processes, and more generally to investigate the molecular basis of host range and susceptibility, virus infectivity and virulence, and the ecology and evolution of baculoviruses.

Figure 1. Assessing assembly quality based on lengths and number of contigs.

(A) Comparison of total numbers of sequences (longer than 400 bps), generated by ABySS or Oases at different k-mer length parameters. The total number of genes from the model insect, B. mori, is shown for comparison. Processed reads were used to run Oases (trimming of 10 last bps and removing of ambiguous quality indicators) and ABySS (no trimming, but removing reads with ambiguous quality indicators) at different k values ranging from 20 to 95 for ABySS and 21 to 85 for Oases. (B) Length distribution of the best Oases and ABySS assembled datasets compared to the B. mori database. Numbers of sequences were classified according to length ranges and plotted for comparison. All EST sequences of B. mori were downloaded from the silkworm database (http://www.silkdb.org).

Figure 2. An example demonstrating assembly coverage by mapping reads to a reference scaffold containing assembled sequences.

The BWA tools mapped original 100 bp reads to the scaffold of all ABySS_55 sequences with criteria that the matches had > = 95% length and > = 90% identity. In this Figure, the top panel shows reference transcripts. The bottom panel shows coverage as the number of reads per nucleotide position in the reference genes. As shown in this example, many parts of the assembled sequences had coverage higher than 200.

Figure 3. BLASTX similarity distribution and top-hit species distribution.

BLASTX was applied for all 29,586 assembled sequences using BLAST2GO (www.blast2go.org/), . (A) Distribution of numbers of sequences at different BLASTX identities. (B) Numbers of top hit sequences from BLASTX were calculated for each species.

Figure 4. Comprehensiveness of H. zea transcriptome as reflected by gene ontology.

Gene ontology terms (GOs) for H. zea were extracted from InterproScan and were compared to those for B. mori downloaded from the silkworm database (http://www.silkdb.org). The GOs were grouped using CateGorizer tool (http://www.animalgenome.org/tools/catego/), .

Figure 5. Relative expression of insect and virus genes from 0–72 hours post infection.

The expression of a gene at a time point is shown in relative level compared to the reference expression of the 28 S RNA gene at that time . Time course expression level for two early H.armigera virus genes, namely DNA polymerase (Pol) and Immediate early gene (IE1) and the host heat shock protein 70 gene (HsP) are shown.

Figure 6. Correlation and distribution of microarray signals.

Correlations of microarray signals between two replicated samples for (A) Infected samples 1 and 2 and (B) un-infected samples 1 and 2 were computed using Limma package in R software, . Each spot represents log 2 signal of the same gene in two samples. The high correlation demonstrates the reliability of detection of the microarray platform used. The distribution patterns of gene ranks based on expression levels of all insect-virus genes (C) or insect genes only (D) are shown for assessing the effect of removing virus genes prior to the normalization step.

Figure 7. Overview of differentially expressed genes.

(A) and (B) show the heatmaps that describe differences in normalized log signal intensity for the first 100 insect genes and 100 virus genes respectively. I1 and I2 are infection 1 and 2, while U1 and U2 are for uninfected samples 1 and 2 respectively. (C) and (D) show overall differential expression profiles for insect genes and virus genes, respectively. The scatter plots show log fold change of expression between two replicates of infected vs. two replicates of uninfected samples (computed by a linear model in LIMMA) and the corresponding Log-odds values. Log-odds is the natural log of the ratio of the probability for the difference being true to the probability of it being not true, i.e. the higher the value the more confidence of difference.

Figure 8. Comparison of numbers of up-regulated genes and down-regulated genes in each GO group, for infected versus non infected samples.

Numbers of up-regulated genes and down-regulated genes were both divided by the total number of genes in the GO group that these genes belong to. The ratio for down-regulated genes was then subtracted from the ratio of up-regulated genes to get a score for each GO group. A positive score suggests the group is up-regulated. Only four groups, namely stress responses, RNA metabolism, enzyme regulator activity and Golgi apparatus were statistically confirmed as being significantly different.

Figure 9. Gene set enrichment analysis (GSEA) of four immune-response pathways in infected H. zea cells.

A GSEA web-based tool was used (http://www.broadinstitute.org/gsea, version 3.7), . The enrichment score (ES) reflects the degree to which a gene set is overrepresented at the top or bottom of an entire ranked list of genes from the microarray data. A positive ES indicates gene set enrichment at the top of the ranked list (more up-regulated); a negative ES indicates gene set enrichment at the bottom of the ranked list (more down-regulated). Each vertical line in the horizontal axis reflects a gene.