Transcriptome Comparison Analysis of Ostrinia furnacalis in Four Developmental Stages

The Asian corn borer, Ostrinia furnacalis, is one of the most destructive pests of maize and causes huge losses in maize yield each year. In order to characterize the different developmental stages, a high-throughput sequencing platform was employed to perform de novo transcriptome assembly and gene expression analysis for the egg, larva, pupa and adult stages. Approximately 185 million reads were obtained, trimmed, and assembled into 42,638 unigenes with an average length of 801.94 bp and an N50 length of 1,152 bp. These unigene sequences were annotated and classified by performing Gene Ontology (GO), Cluster of Orthologous Groups (KOG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional classifications. Comparison of the gene expression profiles of the two transitional stages revealed dramatic differences. Some differentially expressed genes are associated with digestion, cuticularization olfactory recognition and wing formation as well as growth and development. In total, 12 putative insect development-related genes were identified. Real-time quantitative PCR (RT-qPCR) results and sequencing based on relative expression levels of randomly selected genes confirmed these expression patterns. These data represent the most comprehensive transcriptomic resource currently available for O. furnacalis and will facilitate the study of developmental pathways, cuticularization, wing formation and olfactory recognition.

genome. Many insects have been sequenced, including the ACB [12][13][14][15] and its sibling species, the European corn borer 16,17 . The deeper sequencing coverage of Illumina sequencing combined with de novo transcriptome assembly will facilitate the characterization of genes in the absence of a reference genome. We used Illumina HiSeq 2000 to generate a substantial dataset of transcript reads at different stages of O. furnacalis development. The primary aims of this study were to compare gene expression levels in different developmental stages and create a database of molecular information to identify genes related to ACB development.

Materials and Methods
Insect Rearing. O. furnacalis larvae were reared on an artificial diet at the Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing, China. The larvae were kept indoors at a temperature of 27 ± 1 °C with an L:D photoperiod of 16:8 h and 70-80% relative humidity. Eggs, 3 rd stage larvae, pupae and three-day-old moths were collected and immediately placed in liquid nitrogen, then transferred to a − 80 °C freezer until use.

RNA Isolation, cDNA Library Preparation and Illumina Sequencing for Transcriptome
Analysis. Total RNA from eggs, larvae, pupae and moths was extracted separately using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions, and cDNA library construction and Illumina sequencing were subsequently performed at Gene Denovo Biotechnology Co., Ltd, China. Briefly, approximately 20 μ g of pooled total RNA samples from individual O. furnacalis developmental stages were digested with DNase I (Sigma, USA), enriched for mRNA via purification with oligo(dT) magnetic beads and fragmented into fragments of 100-400 bp. First-strand cDNA was synthesized using SuperScript II and random primers. Subsequently, second-strand cDNA was synthesized using DNA polymerase I and RNaseH. The cDNAs were fragmented, end-repaired, and ligated with library-specific barcoded adaptors. Library fragments were PCR-amplified for a minimum number of cycles to avoid normalization for downstream quantitative gene expression analyses. Amplified products were purified with a QIAGEN MiniElute PCR Purification Kit (Qiagen, Venlo, Netherlands), and approximately equal molar proportions of male and female indexed libraries were sequenced in a single flow cell of an Illumina HiSeq 2000 sequencing machine.
Assembly and Functional Gene Annotation. The raw data outputs from the Illumina equipment were trimmed for quality scores (q) < 20, and adaptor, N < 10% and low-quality reads were removed to obtain high-quality, clean reads. The clean reads were assembled using Trinity software (http://trinityrnaseq.sourceforge. net/) to produce unigenes. Subsequently, overall unigenes were analyzed and annotated. Unigenes were aligned with the Nr, Nt, SwissProt, and TrEMBL databases using BLAST with a cut-off E-value of 10 −5 . Functional gene annotations were collected for all unigene sequences ≥ 150 bp using Blast2GO 18 . Initial searches of the National Center for Biotechnology Information (NCBI) non-redundant (Nr) protein database were conducted using the BLASTx algorithm, followed by the collection of Gene Ontology (GO) terms from the GO database and the retrieval of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway designations. Differential Gene Expression and Trend Analysis. Differential unigene abundances were determined by performing independent alignments of short reads from four libraries against the set of O. furnacalis lifetime unigenes using Blat software. Reads per kb per million reads (RPKM) were calculated as: RPKM = (1000000*C)/ (N*L/1000), where C is the number of mappable reads in specific unigenes, N is the total number of reads mapped to unigenes in a particular sample, and L is the length of the unigene 19 . To analyze the expression profiling of eggs, larvae, pupae and moths based on RPKM values, Short Time-series Expression Miner (STEM) software (http://www.cs.cmu.edu/~jernst/stem) was used to compare the trends exhibited in these four stages 20 . P-values correspond to the differential gene expression test, which was performed to analyze all trends in these four stages.
O. furnacalis unigene sequences that mapped to the reference canonical pathways in the KEGG database were analyzed. We specifically assigned 14,947 unigene sequences to 236 KEGG pathways; 1,364 unigenes were involved in metabolic pathways, and these pathways played dominant roles in the following pathways: RNA transport, insect hormone biosynthesis, olfactory transduction, ubiquitin-mediated proteolysis and others (Table S2).

Differential Gene Expression among O. furnacalis Developmental Stages.
To compare the differential expression levels of genes between the major developmental stages, up-and down-regulated gene numbers were calculated between every pair of O. furnacalis life-stages that represented specific transitional developmental stages in holometabolic insects (eggs and larvae, larvae and pupae, pupae and adults, adults and eggs). The top 10 genes for each comparison are listed in Table S3. Based on analyses of the pathways, annotation and the reported literature, 12 insect development-related genes were identified and submitted to the NCBI database.
In total, 19,401 different genes demonstrating significant expression profile changes in the hatching-to-larva process, including 8,168 up-regulated genes and 11,233 down-regulated genes, were obtained in the comparative analysis between eggs and larvae. Four of the ten most up-regulated genes in the egg transcriptome had homology to B. mori genes of neuropeptide-like gene 4, mucin-5AC-like gene, salivary cysteine-rich peptide precursor gene, and cuticular protein gene RR-1 motif 20 isoform X1. Four genes with unknown functions were found to be similar to GJ14248 of Drosophila virilis, CG12017 of Papilio xuthus, uncharacterized protein gene LOC105388107 of Plutella xylostella and hypothetical protein KGM_22714 of Danaus plexippus. Two genes had Scientific RepoRts | 6:35008 | DOI: 10.1038/srep35008 putative functions as a secreted protein of Papilio polytes and a serine protease of Bombyx mandarina. Among the ten most down-regulated genes, three had predicted functions as the proteoglycan, mucin and segmentation proteins of P. xylostella. Two genes putatively matched the pair-rule protein and leucine-rich repeat extensin-like protein 5 of B. mori. Two genes matched the hypothetical protein KGM_13352 and glucose dehydrogenase of D. plexippus. One gene was predicted to function as the epidermal growth factor receptor of Bicyclus anynana and one as the histone 2A variant of Spodoptera frugiperda (Table S3).
In the larva-to-pupation process, comparison between the larval and pupal transcriptomes revealed 20,130 significantly differentially expressed genes, including 6,251 up-regulated genes and 13,879 down-regulated genes in the larvae library. Three of the top ten up-regulated genes were aligned with B. mori, with the predicted functions of larval cuticle protein LCP-22 precursor, zonadhesin-like isoform X2 and uncharacterized protein LOC101735588. Two genes matched with two hemipteran insects (Monomorium pharaonic and Linepithema humile), with functions as collagen alpha-5(IV) and mucin-19. Furthermore, two genes matched with two Noctuidae insects (Spodoptera exigua and Helicoverpa armigera), with MG17 and Niemann-Pick type C2 protein functions. Two genes that matched closely related species (O. nubilalis and Cnaphalocrocis medinalis) were predicted to have the functions of a trypsin-like serine protease and a chemosensory protein. Only one gene was similar to the adhesive plaque matrix protein in Saccoglossus kowalevskii. The top ten hits for down-regulated genes in the larvae library included the following predicted functional genes: the glucose dehydrogenase gene and the flocculation protein gene in B. mori; the uncharacterized protein gene LOC105388107, the mucin-17 gene and the formin-J-like gene in P. xylostella; the CG12017 gene in P. xuthus; the peritrophin type-A domain protein 3 gene in Danaus plexippus; the seminal fluid protein gene CSSFP002 in Chilo suppressalis; and the collagen alpha-2(I) gene in Columba livia (Table S3).
Differentially expressed genes were observed before and after ACB emergence. A total of 18,833 significantly expressed genes (including 8,752 up-regulated genes and 10,081 down-regulated genes) were obtained in the pupae library. Among the top ten gene hits, five genes matched with P. xylostella with functions of the nose resistant to fluoxetine protein 6-like gene, the uncharacterized protein gene LOC105384628, the uncharacterized protein gene LOC105388107, the probable ATP-dependent helicase gene PF08_0048 and the troponin C gene. The predicted functions of the other five top-ten genes included the nuclear-pore anchor-like gene (Bombyx mori), the insect intestinal mucin 4 gene (Danaus plexippus), the retinol-binding protein gene (Papilio xuthus), the chitin binding domain 3 protein gene (Mamestra configurata) and the general odorant-binding protein 3 gene (Cnaphalocrocis medinalis). Among the top ten down-regulated genes in the pupal transcriptome, five genes had predicted functions as the larval cuticle protein LCP-17 precursor gene (B. mori), the GL14840 gene (Drosophila persimilis), the chitin deacetylase 2 gene (Mamestra configurata), the chymotrypsin-like serine protease gene (O. nubilalis) and the pancreatic triacylglycerol lipase-like gene (B. mori). The other genes with predicted functions were the trypsin serine protease gene (O. furnacalis), the lipase gene (H. armigera), the trypsin-like serine protease gene (O. nubilalis), and the mucin-19-like gene (Linepithema humile) (Table S3).
When differential gene expression was assessed between the adult and egg transcriptomes, 24,901 genes exhibiting significantly differential expression were identified. Of these genes, 17,741 were up-regulated and 7,160 were down-regulated in the adult transcriptome (Fig. 5). Five up-regulated genes matched the uncharacterized protein genes KIAA0753 (Erinaceus europaeus), LOC105381897 and LOC105387542 (P. xylostella) and LOC101744818 and LOC101741925 (B. mori), and five up-regulated genes had predicted functions, specifically the repetitive proline-rich cell wall protein 2-like gene (Chrysochloris asiatica), the hypothetical protein gene KGM_11707 (Danaus plexippus), the flexible cuticle protein 12-like gene (P. xylostella), the Kunitz protein 8 gene (Echinococcus granulosus) and the histone H2A sperm-like gene (P. xylostella). Interestingly, 7 down-regulated genes matched to  (Table S3).

Trend Analysis.
To examine the expression profiles of the 27,986 differentially expressed genes (DEGs), the expression data from eggs, larvae, pupae and adults was clustered into 15 profiles by Short Time-series Expression Miner software (STEM), in which 27,986 were clustered into 7 profiles (p-value ≤ 0.05) (Fig. 6), including genes demonstrating down-regulated patterns (Profile 0), comprising 4,233 DEGs. Profile 7, which contained 3,467 DEGs, showed down-regulation from the larvae to the pupae stage and maintained balance during the transitions from eggs to larvae and pupae to adults. Similar to profile 7, profile 8, which included 1,868 DEGs, demonstrated up-regulation during the emergence stage. Profile 1 showed down-regulation from eggs to pupae and then up-regulation during the transition to adults; this profile included 2,217 DEGs. Profile 2 comprised 2,308 DEGs and showed down-regulation during the hatching stage (eggs to larvae), then did not change during the other stages. Profile 9 contained 1,908 DEGs and showed no changes prior to pupation, then exhibited down-regulation Up-regulated unigenes are marked in red, and down-regulated unigenes are marked in green.

Figure 6. DEG expression profiles for the four developmental stages based on P-values.
Colored squares represent significant enrichment, and non-colored squares represent no significant enrichment. "1" represents the egg stage, "2" represents the larval stage, "3" represents the pupal stage, and "4" represents the adult stage.
during the adult stage. Profile 13 contained 1,203 DEGs and showed up-regulation during the larval stage, then was down-regulated during the adult stage.

RT-qPCR Analysis.
To confirm the accuracy and reproducibility of the transcriptome analysis results, RT-qPCR analyses were performed on eight randomly selected genes: Unigene0040326, Unigene0011488, Unigene0032550, Unigene0001704, Unigene0006956, Unigene0012958, Unigene0021965, and Unigene0012689. β -actin served as the reference gene for RT-qPCR normalization. The analyzed results in Fig. 7 supported the differentially expressed unigenes (DEU). Unigene0040326, Unigene0040326, Unigene0006956, Unigene0021965 and Unigene0012689 were primarily expressed in pupae, whereas Unigene0011488 and Unigene0012958 were primarily expressed in eggs. Only Unigene0032550 was primarily expressed in adults. The high confirmation rate of the unigenes indicated the reliability of the transcriptome data.

Discussion
In this study, we compared the transcriptomes of four developmental stages (eggs, larvae, pupae and adults) of the maize pest O. furnacalis. A profound understanding of the molecular mechanisms regulating pest life cycles and development at each stage may aid in the control of this pest by facilitating the development of more sustainable and environmentally-friendly approaches. These results significantly increase the molecular resources available for the study of other insect pests and provide a framework to understand changes at the gene expression level during insect development. We assembled transcriptomes from all four stages and obtained a total of 42,638 unigenes with an average length of 801.94 bp. Approximately 23,494 sequences from this insect were annotated based on four main databases. Approximately half of the total unigenes had not been previously annotated; thus, our Illumina sequencing depth and analysis in this study represent improvements over previous reports 22,23 . Many assembled unigenes were not significantly matched with available databases due to their short sequences or because they represented significantly novel genes. Seven of the top 10 annotated unigenes belong to Lepidopteran, highlighting the conservation of many genes in Lepidopteran. Although we lack the full genome information for O. furnacalis, annotation of the four developmental stages transcriptomes indicated a high proportion of functional genes for this insect.
The ACB exhibits the typical developmental characteristics of holometabolous metamorphosis, and insects develop from eggs to five-stage larvae, then to pupae, and eventually emerge as moths. Moths lay eggs and begin a new life cycle. During the transition from eggs to larvae, 8,168 up-regulated genes and 11,233 down-regulated genes were identified. Among the ten most up-regulated genes, many were related to immunity, digestion and cuticularization, such as the salivary cysteine-rich peptide precursor gene and the cuticular protein gene RR-1 [24][25][26] . The cuticular protein genes have also been identified in the larval stages of Athetis lepigone and other insects, which indicates these protein genes are ubiquitous in insects 23 . Insect cuticles are composed of cuticular protein and chitin, which not only support and maintain the physical structure of the organism but also serves as natural barriers against external adverse impacts 27 . During the transition from egg to larva, cuticle composition and performance rapidly change due to alterations in cuticular proteins; these changes facilitate the performance of cuticle-based structures 28 . Our results provide useful information for further studies investigation cuticular proteins and the cuticle. Furthermore, the segmentation protein was down-regulated in larvae, suggesting this gene functions in cell mitosis. When larvae and pupae were compared, one gene encoding Niemann-pick type C2 protein was identified; this protein is an essential carrier protein for cholesterol and may function in chemical communication from late endosomes and lysosomes to other cellular organelles 29 . The putative peritrophin type-A domain protein 3 gene was also identified in the transcriptome. The peritrophic membrane (PM) is a film-like structure that separates food from the insect midgut tissue. It protects the epithelium against food abrasion and microorganisms and has other functions involved in enzyme compartmentalization 30 . The peritrophin type-A domain protein 3 genes may also carry out additional functions during feeding stages. At the adult stage, the Figure 7. RT-qPCR analysis of eight randomly selected genes to confirm expression patterns indicated by sequencing. Three technical replicates were performed for each of the three biological replicates. The height of each bar chart represents the mean average of sample-specific 2 −ΔΔCt values. "1" represents the egg stage, "2" represents the larval stage, "3" represents the pupal stage, and "4" represents the adult stage.
Scientific RepoRts | 6:35008 | DOI: 10.1038/srep35008 retinol-binding protein gene, an olfactory related gene, a general odorant-binding protein gene and the troponin C gene were identified. These genes are important for moth flight and the ability of moths to seek mates and hosts, particularly the general odorant-binding protein (GOBP) and the troponin C gene. GOBP is small water-soluble protein that transports hydrophobic odorants through the aqueous sensillar lymph and plays an important role in insect chemoreception 31,32 . The troponin C gene may carry out functions for moth flight. Moth flight behavior is dependent on the wings, which are moved by resonant changes in the shape of the thorax produced by indirect flight muscles; muscle contraction is regulated by changes in the Ca 2+ concentration in the fibers. Troponin C, which has a single Ca 2+ -binding site near the C-terminus, regulates these muscles in many insect genera 33 .
Based on the analyzed pathways, annotations and reported literature, 12 insect development-related genes were identified based on the comparison of the transcriptomes at different stages. Many genes participate in insect developmental processes. Various genes, such as segmentation genes, control cellular identity during early pattern formation in Drosophila 34 . Unigene0011488 was annotated as segmentation protein fushi tarazu, which has putative function in development at the blastoderm stage. Furthermore, this gene is transiently expressed in a specific subset of neuronal precursor cells, neurons (such as aCC, pCC, RP1, and RP2), and glia in the developing central nervous system (CNS), suggesting a role in the transformation of neuronal identity 34 . Unigene0012958 was annotated as an epidermal growth factor receptor; Unigene0001331, Unigene0016506, Unigene0001330, Unigene0021965 and Unigene0026727 were annotated as cuticular proteins. These two gene types have previously been reported to function in epidermal growth 35 , cuticle sclerotization and adult molting 36 . The neuropeptide-like 4 gene (Unigene0033676) is expressed uniquely during diapause, suggesting a potential role in initiating and maintaining diapause in the flesh fly 37 . Other genes, such as chemosensory protein (Unigene0032550) and chitin deacetylase 2, have functions related to larval development 38 and developmental timing delays 39 . All identified genes represent candidate genes for future studies.
This study constructed four transcriptome libraries and generated much-needed genetic and genomic resources for the investigation of O. furnacalis. Future analyses of genes related to insect development or death will be useful for pest control. Additionally, the abundance of unigenes and expression data derived from this study will supply information regarding the identification of genes involved in O. furnacalis development.