Transcriptomic Analysis of Eucryptorrhynchus chinensis (Coleoptera: Curculionidae) Using 454 Pyrosequencing Technology.

Eucryptorrhynchus chinensis Olivier (Coleoptera: Curculionidae) is one of the most important pests of Ailanthus altissima; however, so far, no studies on the genome or transcriptome of E. chinensis have been reported. Using the Roche 454 FLX Titanium platform, an RNA pool obtained from E. chinensis eggs, larva, pupae, and adults was sequenced and assembled de novo to achieve maximum diversity of sampled transcripts. We obtained 1,441,137 (∼518 Mb) raw reads with an average length of 360 bp. After trimming, 89% qualified reads were produced and assembled into 35,509 isotigs with an average length of 440 bp, N50 of 1,048 bp, and 111,643 singletons. We generated 87,894 unigenes following a cluster analysis of the isotigs and singletons, and then functionally annotated the unigenes with gene descriptions. We obtained 23,363 GO assignments, and 12,724 unigenes were assigned to KOG. Based on these annotations, 294 biochemical pathways involved in growth, reproduction, and stress or immune responses were predicted. A total of 659,026 single nucleotide variants and 6,112 simple sequence repeats were detected. Our data provide comprehensive information on the sequence and possible functions of E. chinensis transcripts.

Eucryptorrhynchus chinensis Olivier (Coleoptera: Curculionidae) is one of the most important pests of the tree-of-heaven, Ailanthus altissima Swingle, which is distributed throughout China. Eucryptorrhynchus chinensis larvae feed on the root tissues of Ailanthus trees, thereby disrupting nutrient and water translocation and usually killing the trees within 3–5 years of the initial attack. Eucryptorrhynchus chinensis and E. brandti, whose larvae feed on the phloem and cambial, respectively, often occur concurrently. Eucryptorrhynchus chinensis and E. brandti are currently considered as potential biological control agents for A. altissima (Ding et al. 2006a). Recent research has focused on morphology and ecology, and only a few genetic studies have been reported (Herrick et al. 2011, McAvoy et al. 2014).

Second- or third-generation sequencing technology has been employed to generate transcriptome data for non-model species, and to provide valuable genetic in3formation even in absence of genomic sequences (Morozova et al. 2009, Wheat 2010). In multicellular organisms, not all genes are actively transcribed in every cell; therefore, by collecting different types of cells or tissues and sequencing the pooled transcripts, a better understanding of transcriptional activity in those organisms can be achieved. Although the similarities and differences between E. chinensis and E. brandti in terms of morphology, behavior, and ecology allow us to understand the phylogeny, speciation, and evolution of E. chinensis (Ding et al. 2006b, McAvoy et al. 2014), genetic data for E. chinensis permits a better understanding of this species in terms of its genetic evolution and biological functions. However, currently, the genomic sequence resources available for beetles are limited (Park et al. 2008, Pauchet et al. 2009, Van Belleghem et al. 2012).

In this study, we used 454 pyrosequencing and de novo assembly to analyze the transcriptome of E. chinensis. The aim of the study was to achieve maximum diversity of transcripts and produce a resource for large-scale gene discovery, which may facilitate studies of gene expression and function (Ekblom and Galindo 2011, Zhu et al. 2012, Scully et al. 2013). The data derived from this study provide a basis for further investigation of the phylogenetic and speciation processes of E. chinensis and closely related species E. brandti.

Materials and Methods

Insect Sample Preparation and RNA Extraction

To obtain maximum diversity of transcripts, mRNA from different tissues across sexes and developmental stages was isolated and pooled.

Samples for RNA extraction were prepared from eggs, larvae, prepupae, pupae, and adults (male and female) (Supp Table 1 [online only]). Adults were collected from a farmland shelter-forest in Lingwu City (N 38º03′, E 106º20′), Ningxia Hui Autonomous Region in northwest China. Eggs laid by adult females were hatched and bred in our laboratory. During the breeding process, beetle samples were collected at different developmental stages, including eggs, different instars of larvae, prepupae, and pupae. Overwintering larvae (4th and 5th instars) were collected from the outdoor soil in November. All samples were immediately placed in liquid nitrogen and then stored at −80 °C until use.

Total RNA from these frozen individuals was isolated using TRIZOL reagent (Invitrogen, Carlsbad, CA) according to the manufacturer’s instructions. For eggs and young larvae, a mixture of multiple eggs (>100) and nymphs of multiple stages (>25) was pooled to provide samples for RNA extraction. Equal quantities of high-quality RNA from the samples were pooled for cDNA synthesis.

cDNA Library Construction and 454 Sequencing

cDNA libraries were constructed using a SMART PCR cDNA synthesis kit (Clontech, Palo Alto, CA). Normalization of the cDNA was conducted using a Trimmer-Direct cDNA normalization kit (Evrogen, Moscow, Russia). cDNA samples were sheared by sonication to produce fragments of ∼100–800 bp, which is an appropriate fragment size range for 454 sequencing. Finally, all libraries were combined into a single pool.

Assembly and Annotation

After a full 454 GS-FLX run, raw reads were generated and the quality of the reads was assessed using FastaQC (version 0.10.1). The raw reads were pre-processed to obtain cleaned and qualified reads. The process by which the data were prepared for downstream processing applications, such as de novo assembly using Newbler (version 2.6), was designed specifically for assembling sequence data generated by the 454 GS-series of pyrosequencing platforms (Kumar and Blaxter 2010).

Unigene Generation and Annotation

The assembled isotig and singleton sequences were clustered and compared using CD-HIT (version 4.5.6) (Li and Godzik 2006). Sequences with >95% similarity were clustered into one class, and the longest sequence of each of the classes was treated as a unigene during later processing.

To determine the functional annotation of unigenes, the unigene sequences were compared with those deposited in public databases using BLASTx (2.2.26+) with a significance threshold of E-value < 1e-5 and identity of >30% (Camacho et al. 2009). Based on the BLASTx results, unigenes were annotated with Gene Ontology (GO) terms using Blast2GO (Conesa et al. 2005). To analyze unigene-relevant biochemical pathways, KEGG maps and KEGG Orthology classification numbers (KO number) were built for pathway analysis using the online KEGG Automatic Annotation Server (KAAS), http://www.genome.jp/kegg/kaas/ (accessed May, 2015) (Moriya et al. 2007). Unigenes were compared with the KOG library using BLASTx and then mapped to the KOG classification (Tatusov et al. 2003).

Identification of SSRs, SNVs, and Indels

MISA (MIcroSAtellite identification tool) program v1.0 software (http://pgrc.ipk-gatersleben.de/misa/ accessed May, 2015) was used to identify and localize potential SSR markers (Thiel et al. 2003). ssahaSNP software was used to detect SNVs and indels in the assembled isotigs (van Oeveren and Janssen 2009). The parameters used in this study are presented in Additional File 1 (online only).

Results

454 Sequencing and Assembly

After a full 454 GS-FLX run, 1,441,137 (∼518 Mb) raw reads were generated. After quality control assessment, 1,277,554 (88.6%) high-quality reads, which ranged from 45 to 1,467 bp with an average length of 372 bp, were obtained (Table 1 and Supp Fig. 1a [online only]). The sequence data were deposited in the NCBI Sequence Read Archive under accession number SRX719565. Assembly of the E. chinensis transcripts resulted in 35,509 isotigs and 111,643 singletons (Table 1).

Table 1.

Summary of 454 Transcriptome Sequencing and Assembly of E. chinensis

Raw results (after trimming)		Assembly results
Total number of reads	1,277,554	Total number of isotigs	35,509
Total read length (bp)	475,325,044	Total isotig length (bp)	26,938,516
Minimum read lengths (bp)	45	Isotig N50 (bp)	1,048
Median read length (bp)	371	Maximum isotig length (bp)	9,139
Maximum read length (bp)	1,467	Mean depth	32
Mean read length (bp)	372	Number of singletons	111,643
GC content (%)	42.03	Total number of unigenes	87,894

The length of the isotigs ranged from 45 to 9,139 bp with an average length of 758 bp and N50 of 1,048 bp. Most of the isotigs (70.3%) were between 100 to 1,000 bp in length, and 8,527 isotigs had a length >1,000 bp (Table 1 and Supp Fig. 1b [online only]). The remaining 111,643 high-quality reads were retained as singletons and ranged from 45 to 1418 bp with an average length of 282 bp. Among these singletons, 24,529 (22%) were >500 bp.

Exogenous Transcripts

Because samples of extracted RNA were obtained from the whole individual, it is possible that some of the reads were derived from endosymbionts commensal with E. chinensis. The trimmed reads were compared against the NCBI database using Blastn with an E-value of 1e-10 and similarity of >80%, to identify the exogenous transcripts.

Significant BLAST hits were obtained for 111,914 (8.8%) trimmed reads, and best matched 2,294 species. More than 50% of the total reads were matched to 13 species (Supp Fig. 2 [online only]). As expected, most of these were insect species, among which Coleoptera were the most abundant; Dendroctonus ponderosae (26,283; 23.5%), Tribolium castaneum (5,014; 4.5%), and Monochamus alternates (3,706; 3.3%) were the top three species hits. However, Wolbachia sp. wRi (2,089; 1.9%) sequences were also highly represented, indicating that E. chinensis may be infected with Wolbachia.

Functional Annotation of the Transcriptome

Unigenes

A cluster analysis of 35,509 isotigs and 111,643 singlets generated 87,894 unigenes.

For annotation of the unigenes, the unigene sequences were first blast-searched against the public databases using BLASTx. This search matched 24,039 unigenes (27.35% of all unigenes) in the non-redundant databases. The species matched in the non-redundant database were classified mainly into Coleoptera, Hymenoptera, and Lepidoptera. The unigenes of E. chinensis showed significant matches with the sequences of Tribolium castaneum (11,580; 48.3%) followed by Dendroctonus ponderosae (4,840; 20.2%), and Curculio glandium (424; 1.8%) (Supp Fig. 3 [online only]).

Annotation

Based on the results of BLASTx annotation, 23,363 GO assignments were generated and divided into three categories, 61.2% for biological processes (14,299), 15.9% for cellular components (3,719), and 22.9% for molecular functions (5,345), which were further categorized into 64 functional groups (Fig. 1).

Fig. 1.

Gene ontology (GO) distribution for the E. chinensis transcriptome The unigenes were annotated in three main categories: biological process, cellular component, and molecular function.

In total, only 12,724 out of 87,894 (14.5%) unigenes were aligned to the KOG database to predict and classify possible functions. These genes were classified into 25 different functional classes, among which, the cluster for translation, ribosomal structure, and biogenesis (2,499; 17.9%) was the largest group (Fig. 2).

Fig. 2.

KOG function classification of the E. chinensis transcriptome In total, 12,724 unigenes were classified into 25 different functional classes. Note: [S] Function unknown; [Z] Cytoskeleton; [Y] Nuclear structure; [W] Extracellular structures; [V] Defense mechanisms; [U] Intracellular trafficking, secretion, and vesicular transport; [T] Signal transduction mechanisms; [R] General function prediction only; [Q] Secondary metabolite biosynthesis, transport, and catabolism; [P] Inorganic ion transport and metabolism; [O] Post-translational modification, protein turnover, chaperones; [N] Cell motility; [M] Cell/membrane/envelope biogenesis; [L] Replication, recombination, and repair; [K] Transcription; [J] Translation, ribosomal structure, and biogenesis; [I] Lipid transport and metabolism; [H] Coenzyme transport and metabolism; [G] Carbohydrate transport and metabolism; [F] Nucleotide transport and metabolism; [E] Amino acid transport and metabolism; [D] Cell cycle control, cell division, chromosome partitioning; [C] Energy production and conversion; [B] Chromatin structure and dynamics; and [A] RNA processing and modification.

To identify the biological pathways active in E. chinensis, the KEGG database was used to categorize gene functions with an emphasis on the biochemical pathways. Ultimately, 14,788 unigenes were assigned to 3,675 KO assignments and 293 KEGG pathways. The 14,788 unigenes with KEGG annotations were classified into five main clusters, including metabolism, genetic information processing, organismal systems, cellular processes, and environmental information processing (Fig. 3).

Fig. 3.

KEGG pathway classification of the E. chinensis transcriptome In total, 14,788 unigenes were classified into five pathways including Metabolism, Genetic information processing, Organismal systems, Cellular processes, and Environmental information processing.

Using the MISA program, a search for microsatellites in the assembly reads revealed 6,112 SSRs in 6,009 sequences. Among the various classes of SSRs, dinucleotide (473) and trinucleotide (957) motifs, which are used frequently as molecular markers in the field of molecular evolution and phylogenetics, were found and the most common repeat motifs were (AG)n and (AAG)n.

More than 659,026 SNVs and indels were identified (488,872 SNVs and 170,154 indels) from the isotigs using the ssahaSNP program (Supp Fig. 4 [online only]). The predicted SNVs included 281,444 transitions and 207,428 transversions. This corresponds to a transition: transversion ratio of 1.36:1. The most common substitution was A (T) to G (C), representing 57.6% of all substitutions. The overall density of all SNVs and indels in the transcriptome was 1 per 89 bp.

Discussion

We used RNA samples isolated from different developmental stages to characterize the transcriptome of E. chinensis using 454 pyrosequencing technology and de novo assembly. A mixed cDNA library was produced by pooling the RNAs to provide a broad and representative pool of transcripts (Wheat 2010).

Comparing our results with those from recently published studies, revealed that the 454 sequencing depth and assembly efficiency in this study were appropriate (Hale et al. 2009, Meyer et al. 2009, Scully et al. 2013, Wang et al. 2013). In this study, the proportion of sequences without annotation information in the published databases was considerable, at ∼73%. This poor annotation efficiency could be due to the lack of phylogenetically related species currently sequenced with data in public databases. On the other hand, for those species whose genomes or transcriptomes have been sequenced and deposited in the databases, annotation efficiency was higher than those without. For example, the common unigenes of E. chinensis showed significant matches with Tribolium castaneum, although these organisms exhibit many differences in their biological and ecological traits. The results indicated that many unigenes may be novel sequences and a high number of Coleoptera- or species-specific transcripts. In contrast, the efficiency of the BLAST comparison depends, in part, on the length of the query sequence, as the length of the sequences increases, the annotation rate gradually increases. For those sequences >1 kb, the proportion of annotation increased to 50%.

Exogenous Genes

More than 80% of the aligned sequences had a significant BLAST hit that best matched insect species, of which Coleoptera were the most abundant. This finding suggests that the vast majority of our sequences originated from E. chinensis, although some transcripts were exogenous.

Investigation of E. chinensis xenobiotics could thus provide information regarding coevolution, disease, and infection. Interestingly, many reads matched Wolbachia, which are symbiotic bacteria known to affect population dynamics and are hypothesized to be present in E. chinensis populations (Charlat et al. 2007).

Wood-feeding insects often work in collaboration with microbial symbionts to degrade lignin biopolymers and release glucose and other fermentable sugars from carbohydrates, including cellulose and hemicelluloses of the recalcitrant plant cell wall (Scully et al. 2013).

Marker Identification and Characterization

The transcriptomic data obtained by 454 sequencing provided an excellent source for mining and developing gene-associated markers.

The SNV frequency in our assembly reads (1.1/100 bp) was found to be higher than that reported in other studies using cDNA pooled from multiple individuals (Morin 2004, Wondji et al. 2007). Compared with equivalent substitution, a biased 1.36:1 transition (ts) to transversion (tv) ratio was identified in our study

In the present study, SNVs were identified from 4,000 isotigs that had a coverage of 10 or more reads, suggesting that these SNVs were covered at a sufficient sequencing depth and were more likely to represent true SNVs. Among the SNVs, 176,330 (36%) were identified from isotigs with GO annotation information. These SNVs are priority candidates for marker development and should be useful in further studies of important traits, such as phenotypes of adaptive traits (Luikart et al. 2003).

In total, 6,116 SSRs and 659,026 SNVs and indels were identified from the assembled sequences and the total abundance and frequency of the SSRs were 4.2% and 0.1/kb, respectively. Of 6,009 SSR-containing sequences, 1,430 (23.8%) sequences had been annotated. Approximately 5,131 (84%) SSR motifs had >50 bases on both sides.

Using the SSR Locator (Carlos da Maia et al. 2008), 20 primers were designed based on the unigenes including SSR at random, and verified by PCR analysis for E. chinensis and E. brandti. After PCR and electrophoresis, 17 (85%) markers amplified DNA fragments in E. chinensis and 13 (65%) in E. brandti (Supp Fig. 5 [online only]). It is feasible to develop SSR markers from unigenes generated by transcriptome sequencing. These markers provide an effective method of analyzing gene polymorphisms and function of E. chinensis and E. brandti in the further studies.

Candidate Gene Identification

Eucryptorrhynchus chinensis is one of the most problematic pests of Ailanthus in China, and is a potential biological control agent for Ailanthus in North America. Although the biology and ecological traits have been well studied, the molecular basis underlying these traits remains unclear.

GO, KOG, and KEGG analyses are valuable tools that can be used to study the molecular basis of including growth, reproduction, stress, and immunity in E. chinensis.

Transcripts of genes putatively involved in growth (GO: 0040007) and reproduction (GO: 0000003) were found in our 454 data sets. GO terms related to olfaction were also present, such as “odorant binding (GO: 0005549)”, “signal transducer activity (GO: 0004871)”, and “response to stimulus (GO: 0050896)”. Those transcripts can be used to examine the co-evolution and interactions between insects and their environment (Hall et al. 2013).

Unigenes with enzymatic activity, such as “hydrolase activity” and “transferase activity,” were well represented. Ryanodine-sensitive calcium-release channel activity (GO: 0005219) is related to insecticide resistance. Heat shock protein (HSP) binding (GO: 0031072) functions to protect cells from the effects of heat and other stress factors. Carboxylesterase activity (GO: 0004091) is hypothesized to play a pivotal role in detoxifying host tree defense chemicals (Li et al. 2007, Zimmer et al. 2014).

Wood-feeding beetles, especially in the larval stage, feed primarily on the heartwood. These beetles must overcome challenges of digesting intractable woody tissue to acquire sufficient nutrients to complete development. In the present study, unigenes with GO annotations involved in carbohydrate metabolism (GO: 0044262, GO: 0005975) were detected. These unigenes will facilitate the identification of transcripts encoding digestive enzymes or microbial-derived enzymes (Pauchet et al. 2009, Scully et al. 2013).

In this study of the transcriptome of E. chinensis, we obtained abundant transcripts related to molecular mechanisms. To further understand the adaptability of E. chinensis, we analyzed the transcripts related to detoxification, Digestion and stressful conditions. we obtained abundant transcripts related to detoxification mechanisms (Additional File 2 [online only]). Three major multigene enzyme families primarily responsible for xenobiotic metabolism were identified in the transcriptome data. We identified 111 unigenes encoding P450s, which were the most abundant sequences among the detoxification-related genes. Furthermore, abundant sequences encoding carboxylesterases (14) and GSTs (31) were also identified. For P450s, the number of unigenes belonging to the CYP4 clade was higher than that in CYP6 and CYP9 clades, which was similar to the number in other insects, in which most P450 genes were grouped into the CYP3 and CYP4 clade (Meng et al. 2015). These findings will be useful for studies on insecticide action, selectivity, and detoxification in E. chinensis and E. brandti. Eucryptorrhynchus chinensis and E. brandti, whose larvae feed on the cambial and root tissues, respectively, often occur concurrently. Digestion-related genes were identified in the E. chinensis transcriptome dataset (Additional File 3 [online only]), particularly, those encoding enzymes capable of degrading cellulose and hemicelluloses, which are the two most predominant polysaccharides found in hardwoods (Scully et al. 2013). Five transcripts of xylanase and 34 of glucanase were detected in the transcriptome assembly, which may derive from endosymbionts commensal with E. chinensis. This finding suggests that collaboration with microbial enzymes may be required for survival in woody tissue. Heat shock protein reflects the adaptive capability of the E. chinensis to diverse environment. In the present transcriptome, we obtained a number of unigenes encoding HSPs (Additional File 4 [online only]) and speculated that these may help E. chinensis to adapt to the varied environments and biological stresses. We identified 40 hits encoding HSP70, and 48 hits encoding HSC70. The heat shock cognate (HSC) protein binds to nascent polypeptides to facilitate correct protein folding. This transcriptome analysis generated a large number of genes newly identified in E. chinensis, representing a substantial contribution to existing sequence resources. Further studies are needed to validate the functions and expression patterns of these candidate genes, and to investigate their potential roles in adaptability to the environment.

Ethics Statement

Eucryptorrhynchus chinensis Olivier (Coleoptera: Curculionidae) is not an endangered or protected species and has recently become an important pest of A. altissima in China. No specific permits were required for collection in these locations or for these activities.

Supplementary data

Supplementary data are available at Journal of Insect Science online.