De Novo Transcriptome Analysis of Wing Development-Related Signaling Pathways in Locusta migratoria Manilensis and Ostrinia furnacalis (Guenée)

Background Orthopteran migratory locust, Locusta migratoria, and lepidopteran Asian corn borer, Ostrinia furnacalis, are two types of insects undergoing incomplete and complete metamorphosis, respectively. Identification of candidate genes regulating wing development in these two insects would provide insights into the further study about the molecular mechanisms controlling metamorphosis development. We have sequenced the transcriptome of O. furnacalis larvae previously. Here we sequenced and characterized the transcriptome of L. migratoria wing discs with special emphasis on wing development-related signaling pathways. Methodology/Principal Findings Illumina Hiseq2000 was used to sequence 8.38 Gb of the transcriptome from dissected nymphal wing discs. De novo assembly generated 91,907 unigenes with mean length of 610 nt. All unigenes were searched against five databases including Nt, Nr, Swiss-Prot, COG, and KEGG for annotations using blastn or blastx algorithm with an cut-off E-value of 10−5. A total of 23,359 (25.4%) unigenes have homologs within at least one database. Based on sequence similarity to homologs known to regulate Drosophila melanogaster wing development, we identified 50 and 46 potential wing development-related unigenes from L. migratoria and O. furnacalis transcriptome, respectively. The identified unigenes encode putative orthologs for nearly all components of the Hedgehog (Hh), Decapentaplegic (Dpp), Notch (N), and Wingless (Wg) signaling pathways, which are essential for growth and pattern formation during wing development. We investigated the expression profiles of the component genes involved in these signaling pathways in forewings and hind wings of L. migratoria and O. furnacalis. The results revealed the tested genes had different expression patterns in two insects. Conclusions/Significance This study provides the comprehensive sequence resource of the wing development-related signaling pathways of L. migratoria. The obtained data gives an insight into better understanding the molecular mechanisms involved in the wing development in L. migratoria and O. furnacalis, two insect species with different metamorphosis types.


Introduction
Insects are the only group of invertebrates that have evolved flight [1]. Their wings serve not only as organs of flight, but also may be adapted variously as protective covers [2], thermal collectors [3], gyroscopic stabilizers [4], sound producers [5], or visual cues for species recognition and sexual contact [6]. For those insects with wings during the adult stage, complete wings are not always visible throughout the life cycle. The insects undergoing incomplete metamorphosis have represented functional wings during the stage of nymph [7], while the insects going through complete metamorphosis only have wing discs inside the body in the larval stage [7]. Therefore, the comparative studies on how the insect wings are developed will be helpful to understand the insect metamorphosis. Additionally, given the high evolutionary conservation of the proteins involved in the wing development, the understanding of the molecular mechanisms involved in wing development also sheds light on organogenesis, tissue homeostasis, human disease and so on [8][9][10][11].
Insect wing development is controlled with amazing precision and complication. Current understanding of insect wing development mechanisms is mainly from the fruit fly Drosophila melanogaster [12][13][14]. The adult wing of fruit fly is derived from the wing imaginal disc, formed at the end of embryonic development [15]. In Drosophila, the wing imaginal disc is subdivided into anterior (A) and posterior (P) compartments at very early stage and then further subdivided into dorsal (D) and ventral (V) compartments at second instar stage [16][17][18]. Organizers located in the A/P and D/V boundaries coordinate the patterning of the wing disc by secreting signal molecules including the long-range morphogens Decapentaplegic (Dpp; the vertebrate homolog of which is TGFb) and Wingless (Wg; the vertebrate homolog of which is Wnt) [11,19,20], and short-range morphogen Hedgehog (Hh) [21]. The produced morphogens form gradients to regulate the expression of target genes and control all aspects of wing development at cell level via the specific signaling pathways [19,[22][23][24][25][26][27]. These signaling pathways are highly evolutionarily conserved in many different animals. Some key genes for the wing disc development have also been identified in a few limited insect species, such as Tribolium castaneum [28]. However, knowledge about the identification and involvement of wing development-related genes in various insects, especially in non-model insects, is still unclear and incomplete.
Orthopteran migratory locust, Locusta migratoria manilensis, and lepidopteran Asian corn borer, Ostrinia furnacalis (Guenée), are two different types of insects undergoing incomplete and complete metamorphosis, respectively [29,30]. Identification of candidate genes regulating wing development would provide insights into the further study about the molecular mechanisms controlling metamorphosis development. Traditionally, such gene identification in ''non-model'' insects relied on degenerate PCR, which is labor-intensive, expensive, prone to failure, and only produces incomplete fragments [31]. The introduction of novel high throughput sequencing technologies greatly facilitates the global analysis of the blueprint of development-related genes [32]. This technology has been used, for example, to characterize the wing development-related genes in the milkweed bug Oncopeltus fasciatus [33], the oriental fruit fly Bactrocera dorsalis [34], and the salt marsh beetle Pogonus chalceus [35] etc.
In this study, we combined the Illumina sequencing and de novo assembly to obtain and characterize the transcriptome of the wing discs of L. migratoria nymph. 91,907 unigenes were assembled and 23,359 ones were annotated to known databases. We also recharacterized the previous transcriptome of O. furnacalis larvae [36], with special emphasis on wing development-related genes. Overall, we identified 50 potential wing development-related unigenes from L. migratoria transcriptome, and 46 unigenes from O. furnacalis transcriptome, respectively. Additionally, we performed qRT-PCR analysis to investigate the gene expression profiles of several key wing development genes during the stage of rapid growth in L. migratoria and O. furnacalis. All these results provide valuable information for studying the molecular mechanisms involved in the insect wing development, and are useful resources for further exploring the mechanism how signaling pathways control wing development.

Dissection and total RNA extraction of wing discs
The wing discs of L. migratoria fifth instar nymph were dissected along the wing root with small scissors under microscope. Thirty pairs of dissected wing discs were combined, and total RNA samples were prepared using TRizol Reagent (TIANGEN, Beijing, China) following the manufacturer's instructions. Total RNA was dissolved in H 2 O and RNA quantity was determined on a Nanodrop ND-2000 spectrophotometer (NanoDrop products, Wilmington, DE, USA). RNA integrity was checked on Agilent 2100 BioAnalyzer (Agilent Technologies, Englewood, CO, USA). Library construction and Illumina sequencing Ten mg of total RNA was used to isolate mRNA using oligo(dT) magnetic beads. The cDNA library was constructed using NEBNext mRNA Library Prep Reagent Set (NEB, Ipswich, MA, USA) following the manufacturer's protocols. Briefly, enriched poly(A) RNA of each sample was fragmented into 200-700 nt pieces with RNA Fragmentation Reagents. The cleaved RNA fragments were transcribed into the first-strand cDNA using random hexamer-primers, followed by second-strand cDNA synthesis. The resulting double-stranded cDNA (dsDNA) was purified with QiaQuick PCR extraction kit (Qiagen, Hilden, Germany) and dissolved in EB buffer. The purified dsDNA was treated with T4 DNA Polymerase and T4 Polynucleotide Kinase for end-repairing and dA-tailing. After that, they were ligated to sequencing adaptors with barcode using T4 DNA ligase. Finally, fragments with around 200 bp-length were purified with Qia-Quick GelPurify Kit (Qiagen, Hilden, Germany), and used as templates for PCR amplification to create the cDNA library. The library was sequenced on Illumina HiSeq 2000 (Illumina, San Diego, CA, USA) in Baimaike company (Beijing, China).

Assembly and annotation of transcriptomes
Raw reads were filtered to remove low quality reads with Q20 less than 20 and the sequence reads containing adapters and poly-A/T tails. The resulting clean reads were assembled to produce unigenes using the short reads assembling program -Trinity [37]. For functional annotations, we first searched all unigene sequences against various protein databases such as Nr, Swiss-Prot, COG, and KEGG using BLASTX, and then searched nucleotide database Nt using BLASTN, with an E-value cut-off of 10 25 [38]. The BLAST results were used to extract coding region sequences (CDS) from the unigene sequences, and translate them into peptide sequences. When a unigene happened to have no BLAST hits, ESTScan software [39] would be used to determine the sequence direction. In addition, we performed the Gene Ontology (GO) annotations for each unigene with Blast2GO program according to the GO association done by a BLASTX against the Nr database [40,41].
Identification and sequence analysis of wing development-related genes from L. migratoria and O. furnacalis transcriptome The available wing development-related gene sequences from Drosophila were used as references to screen L. migratoria transcriptome database obtained above and O. furnacalis transcriptome obtained previously [36]. The potential candidates of L. migratoria and O. furnacalis wing development-related genes were confirmed by searching the BLASTX algorithm against the non-redundant (nr) NCBI nucleotide database using a cut-off Evalue of 10 25 .
For the sequence analysis of putative wing development-related genes identified above, the deduced protein domains were determined by using Pfam (http://www.sanger.ac.uk/Software/ Pfam/) and SMART (http://smart.embl.de/). Analysis of deduced amino acid sequences, including prediction of signal peptide, molecular weight and isoelectric point, was carried out in the EXPASY (Expert Protein Analysis System) proteomics server (http://www.expasy.org). Sequence comparisons and phylogenetic analysis were performed by MEGA5 software [42]. Phylogenetic trees were constructed by the neighbor-joining

Expression assay of several identified wing developmentrelated genes
To investigate the expression profiles of several key wing development-related genes in forewings and hind wings, we dissected the wing discs as described above, and extracted total RNA independently from 3 biological replicates. DNase I-treated RNA (1mg) was converted into first-strand cDNA using TIAN-Script RT Kit (TIANGEN, Beijing, China). The cDNA products were diluted 2 fold for use as template. Specific primers for each gene were designed and listed in Table S1. L. migratoria Actin and O. furnacalis ribosomal protein L8 (rpL8) was used as an internal standard to adjust the template amounts in a preliminary PCR experiment. The qRT-PCR was performed on a Applied Biosystems 7500 Real-time PCR system (Life Technologies, Grand Island, NY, USA) using the GoTaq qPCR Master (Promega, Madison, WI, USA), according to the manufacture's instructions. The thermal cycling conditions for qRT-PCR and calculation methods were same as described previously [36].

Results and Discussion
Dissection and observation of wing discs in L. migratoria and O. furnacalis All insects in the Pterygota undergo metamorphosis from immature to adult. Among them, insects with incomplete metamorphosis, such as migratory locust L. migratoria, have young nymph resembling the adult with visible forewings and hind wings. Meanwhile, insects with complete metamorphosis, such as corn borer O. furnacalis, go through four stage processes from egg, larva, pupa, and adult in which only adult has visible wings. Although there are some minor differences, the overall appearances of forewings and hind wings in L. migratoria nymph and adult are similar (Fig. 1). However, actual wings are only visible in O. furnacalis adult. In fifth-instar O. furnacalis larva, a pair of pea-like wing discs composed of large amount of cells was dissected from the position where the forewing and hind wing will be derived. There is no similarity on the appearance between the wing discs and adult wings (Fig. 1). The huge difference about the wing development in L. migratoria and O. furnacalis suggests that the molecular mechanisms controlling the wing development in insects with incomplete or complete metamorphosis might be largely different. As a first step to understand these molecular mechanisms, it is important and necessary to identify as many as possible genes functioning in the wing development. Previously, we have combined the Illumina sequencing and de novo assembly to obtain the high quality transcriptome from O. furnacalis larvae [36]. In this study, we also obtained the data for L. migratoria transcriptome, and we can comparatively characterize both transcriptomes and identify wing development-related genes. This work will provide useful information for studying the molecular basis involved in the wing development in L. migratoria and O. furnacalis, two insect species with different type of metamorphosis.

Sequencing and unigene assembly of L. migratoria transcriptome
In order to obtain detailed information about L. migratoria transcriptome, we prepared cDNA from the wing discs of fifthinstar nymph, and subjected it to Hiseq 2000 sequencing. After cleaning of dirty reads and quality checks, a total of 83,000,540 high-quality clean reads (SRA accession number SRX491784) with a cumulative length of 8,382,624,748 nucleotides (8.38 Gb) were generated from L. migratoria wing disc library. The GC percentage of the reads is 42.65%, which is comparable with genome sequence of other insects. Using Trinity software by the manner of paired-end joining and gap-filling, these reads were assembled into 91,907 unigenes longer than 200 nt (mean length of 610 nt and N50 of 1024 nt). The size distribution indicated that the ratio of unigenes with a length of 200-1000 bp was 86.49%, while the length of 12,420 (13.51%) unigenes was more then 1000 bp (Fig.S1). It was significantly larger than that in previous insect transcriptome projects [43,44]. The assembled sequences have been deposited in the NCBI Transcriptome Shotgun Assembly (TSA) Database under the accession GBDZ00000000.

Functional annotation and classification
In order to annotate the unigenes, the obtained sequences were first aligned by BLASTX to various protein databases of nr, Swiss-Prot, KEGG, COG, and GO (E-value,10 25 ), and then aligned by BLASTN to nucleotide database nt (E-value,10 25 ). Among 91,907 unigenes, 20,746 (22.6%), 8,737 (9.5%), 11,450 (12.5%), 5,312 (5.8%), 4,986 (5.4%), and 10,644 (11.6%) ones were annotated in nr, nt, Swiss-Prot, KEGG, COG, and GO, respectively (Table S2). A total of 23,359 unigenes (25.4%) were annotated to at least one database. The remaining 68,548 unigenes (74.6%) were not annotated to any referred databases. The low annotated percentage might be due to the transcripts derived from the cDNA of untranslated regions, assemblage errors, and nonconserved areas of proteins where homology is not detected [43]. Another possibility was that a large part of the genes in L. migratoria transcriptome database were with unknown functions, and the un-annotated unigenes were the potential sources of novel genes.
Gene Ontology (GO) assignment programs were used to classify the functions of all unigenes. GO contains three categories: biological process, cellular component, and molecular function.
Using the Blast2GO and WEGO software, 43,415, 24,188, and 13,706 unigenes were associated with biological process (22 subcategories), cellular component (16 sub-categories), and molecular function (17 sub-categories), respectively (Fig. 3). Among the biological process assignments, cellular processes (15.0%) and metabolic processes (14.4%) represented the most abundant subcategories. It indicated the importance of cell cycle, generation as well as metabolic activities in wing development stage. Under the category of cellular component, the top 3 sub-categories were cell part (21.0%), cell (20.2%) and organelle (16.0%). In the molecular function category, binding (43.0%) and catalytic activities (36.9%) were the most abundant (Fig. 3). The subcategories taking up the largest two proportions in each category were consistent with that in transcriptomic studies of other insects [45,46].
The unigenes were further compared against COG for the analysis of the putative protein functions. A total of 6,508 unigenes were functionally classified into 25 COG categories (Fig. 4). The cluster for 'General function prediction only' (1,407 unigenes) constituted the largest group, followed by 'Replication, recombination and repair' (809 unigenes), 'Translation, ribosomal structure and biogenesis' (560 unigenes), and 'Transcription' (441 unigenes) (Fig. 4). The group of 'signal transduction mechanisms' contains 354 unigenes, the fifth largest one within the COG classification of the L. migratoria transcriptome. It suggests the possible importance of signaling transduction pathways in wing development in L. migratoria.       [20,22,25,[47][48][49]. In order to obtain the global perspective on the wing development-related transcripts in L. migratoria, we searched the assembled transcriptome for orthologous genes known to be involved in the above four pathways in Drosophila.
In total, we have identified 50 unigenes with significant similarity to Drosophila wing development-related genes, including 7 ones in Notch pathway, 12 in Hh pathway, 9 in Dpp pathway, and 10 in Wg pathway (Table 1). Meanwhile we performed the similar analysis for O. furnacalis transcriptome reported previously [36], and identified 46 unigenes with high potential to be related to the wing development (Table 1). Genes involved in the Hh signaling pathway. The Hh pathway is an important signaling cascade in Drosophila to pattern the embryonic cuticles and adult appendages [50]. It is also vital for diverse aspects of animal development and essential in humans to pattern their limbs and internal organs [51,52]. The pathway takes its name from its polypeptide ligand, an intercellular signaling molecule called as Hh [53]. Hh was initially discovered as a segment polarity gene in Drosophila [54]. Mammals have three Hh orthologues, Sonic (Shh), Indian (Ihh), and Desert (Dhh) [55]. We have identified the potential hh genes, a37032 from L. migratoria and Unigene16362 from O. furnacalis transcriptome, respectively (Table 1). Although these two unigenes were incomplete, they have 66% and 73% amino acid sequence identity to Drosophila Hh (Fig. 5). On the other hand, Drosophila Hh is synthesized as inactive precursors with an N-terminal signaling region linked to a C-terminal autoprocessing region which begins with a cysteine for an acyl rearrangement analogous to step 1 of the protein splicing pathway [56]. This Cys residue is completely conserved in both L. migratoria and O. furnacalis Hh homologues (Fig. 5). The conserved motif A, B, F, J, K, and L in Drosophila Hh [57] were found in L. migratoria, but not in O. furnacalis Hh because its C-terminus was presently unknown (Fig. 5). The phylogenetic analysis of Hh from different insect species and other animals reveals that almost all nematode Hh-related proteins form a distinct clade while the other eukaryote Hh proteins, including L. migratoria a37032 and O. furnacalis Unigene16362, cluster into another Hh clade (Fig. 6). L. migratoria a37032 is grouped with other Orthopteran Hhs, and O. furnacalis Unigene16362 is clustered with other Lepidopteran Hhs. We performed qRT-PCR analysis to analyze the transcriptional expression of hh during different developmental stages. The mRNA abundance of hh in L. migratoria remained unchanged from the first instar through the fifth instar stage. The transcript abundance of O. furnacalis hh was consistent before the fourth instar stage, but increased significantly in the fifth instar larvae (Fig. 7). This result also suggests that the quality of the transcriptome assembly is high enough to be used for primer design.
Hh signaling is mediated by a multi-component receptor complex in the cell membrane. This receptor complex consists of a 12-span transmembrane protein, Patched (Ptc) as the receptor and a 7-span transmembrane protein, Smoothened (Smo) as the obligatory signal transducer across the plasma membrane [51,58]. When extracellular Hh binds to and is inhibited by Ptc, Smo starts to accumulate, and inhibit the proteolytic cleavage of zinc-finger transcription factor Cubitus interruptus (Ci) which is normally bound by the kinesin-like protein Costal-2 (Cos2). The intact Ci protein then translocate into the nucleus, allowing the transcription of some genes such as dpp. In the absence of Hh, Ptc blocks Smo activity, and full-length Ci protein is degraded into Ci fragment (CiR) that functions as a transcriptional repressor to block the transcription of target genes (Fig. S3A) [50]. By searching L. migratoria and O. furnacalis transcriptome using Drosophila corresponding genes as reference, we have identified Genes involved in Dpp signaling pathway. Dpp, directed by the Hh signaling, is a member of the bone morphogenetic protein (BMP) growth factor family [59]. Dpp forms a long-range dynamic and precise gradients to control cell survival, cell morphogenesis, cell proliferation, cell differentiation during the wing development [22,25,26,60]. When a cell receives a Dpp signal, its heteromeric receptor complex composed of type I receptor Thickveins (Tkv) and type II receptor Punt is activated [61,62]. Once activated, the receptors are able to phosphorylate an intracellular protein called mothers against Dpp (Mad) [63]. The phosphorylated Mad then associate with the Medea (Med), and the complex translocates to the nucleus where it binds to DNA and activates or suppresses the expression of the target genes in conjunction with other transcription factors [64]. Within this signaling pathway, the Drosophila Hox gene Ultrabithorax (Ubx) restricts both the transcription and the mobility of Dpp [65]. Genes activated by Dpp signaling pathway include transcription factor optomotor-blind (omb) and two members of Spalt (sal) family spalt major (salm), spalt related (salr) and so on (Fig. S3B) [66,67].
In this study, we identified almost all transcripts of the above key components in Dpp signaling pathway (Table 1). Among them, one unigene with sequence similarity to Drosophila Dpp was identified from the transcriptome of L. migratoria (a136595) and O. furnacalis (Unigene7326), respectively. The 59-end of the cDNA sequences of both potential Dpp is incomplete. In order to establish the orthology of these two fragments, we performed a phylogenetic analysis incorporating a selection of BMP growth factor members from Drosophila, other arthropods and a variety of deuterostome taxa. The resulting phylogenetic tree distinguishes two groups of proteins. One group consists of the arthropod dpp genes and their deuterostome homologs, the BMP2/4 genes. The second group comprises the Drosophila TGF-beta genes screw (scw) and glass bottom boat (gbb), and the remaining deuterostome BMP genes with TGF-beta homology (Fig. 8). The identified fragments, L. migratoria a136595 and O. furnacalis Unigene7326, reside in the dpp group with statistical support (reliability value = 98 and 99 for a136595 and Unigene7326, respectively). Therefore, they are postulated as potential dpp genes. We conducted qRT-PCR analysis to investigate the expression profile of dpp. L. migratoria dpp was expressed at high level in the first instar nymphs, then decreased significantly in the second, third, and fourth instar nymphs, but significantly increased in the hind wings of the fifth instar nymph. The mRNA level of O. furnacalis dpp remained unchanged in all tested stages (Fig. 7). The different expression patterns of dpp in two insects suggest that dpp might be related to the different metamorphosis.
In addition, it is interesting that only one mad gene was present in Drosophila but two transcripts encoding for potential mad were identified from both L. migratoria and O. furnacalis transcriptomes (Table 1). Two Mads in L. migratoria (a11173 and a12619) shared 60% identity while two O. furnacalis Mads (Unigene16360 and Unigene6852) had 61% identity in amino acid sequences. A similar situation happened to the identification of med. There is only one med gene in Drosophila but two transcript fragments were identified for both L. migratoria and O. furnacalis med, including a26380 and a11651 in L. migratoria, and unigene4133 and unigene17363 in O. furnacalis, respectively. The encoded amino acid sequences of a26380 and unigene4133 were highly similar to the N-terminus of Drosophila Med while a11651 and uni-gene17363 were highly similar to the C-terminus of Drosophila Med (Fig. S2). We doubted that the two unigenes from each transcriptome were just two partial fragments within med gene. However, no overlapped sequences were observed in a26380 and a11651, or in unigene4133 and unigene17363 (Fig. S2). One possible reason might be due to the sequencing errors because the predicted overlapping part was just located at the terminus of sequenced fragment. Further experiments are required to determine the full length of med in these two insect species.
Genes involved in the Notch signaling pathway. The Notch signaling pathway is highly conserved throughout the animal kingdom. It regulates cell-fate determination during development and maintains adult tissue homeostasis [47]. The key component of this signaling pathway is the Notch receptor. Notch is a 300 kDa single-pass transmembrane protein, composed of a large extracellular domain, a single transmembrane portion, and a small intracellular region [68][69][70]. After binding to its ligands, Delta and Serrate (known as Jagged in mammals), inactive Notch precursor undergoes two proteolytic cleavage events: the first cleavage is catalyzed by ADAM-family metalloproteases; the second cleavage is mediated by c-secretase which is an enzyme complex containing presenilin, nicastrin, PEN2 and APH1 [71]. The second cleavage liberates the Notch intracellular domain (NICD), which then migrates into the nucleus and cooperates with the DNA-binding protein CSL (also named as CBF1, Suppressor of hairless (Su(H)), LAG-1, or RBP) and its co-activator Mastermind (Mam) to promote the transcription of downstream target genes, such as wg (Fig. S3C) [72]. Using the corresponding components in Drosophila as referred sequences, we have identified 9 transcripts for potential components in Notch signaling pathway from L. migratoria and O. furnacalis transcriptomes, respectively (Table 1). Four out of nine unigenes in L. migratoria are complete while only one unigene is complete in O. furnacalis. The possible reason was that the sequencing quality of L. migratoria transcription was better than that of O. furnacalis transcriptome (8.38 Gb vs. 4.72 Gb). Given the importance of Notch in this signaling pathway, we selected it for further analysis. Similar to the case in other invertebrates, we only identified one transcript for Notch gene, a1269 from L. migratoria transcriptome, and Unigene6905 from O. furnacalis transcriptome ( Table 1). The canonical Notch protein consists of three repeated sequence motifs: the extracellular domains contain 10,36 copies of an ,40-amino acid epidermal growth factor-like (EGFL) sequence motif and 3 copies of an ,40-amino acid lin/Notch/glp (LNG) sequence motif; the intracellular domains contain 6-7 copies of a cdc10/SW16/ankyrin (CDC/ANK) sequence motif flanked by stretches of nonrepetitive sequences [73]. L. migratoria a1269 encodes a 2,484-amino acids full-length protein which has 63% similarity to Drosophila Notch. L. migratoria Notch contains 36 EGFL tandem repeats, 3 LNG repeats, and 7 tandem ANK repeats (Fig. 9A), suggesting it is a canonical Notch. O. furnacalis Unigene6905 encodes a 780-amino acids polypeptide with 70% similarity to Drosophila Notch. Only 20 tandem EGFL repeats are predicted in the current identified Unigene6905 fragment (Fig. 9A). It is unknown whether O. furnacalis Notch also contain classic LNG and ankyrin repeats because the current transcript is incomplete. We performed the phylogenetic analysis for L. migratoria and O. furnacalis Notch and 30 Notch sequences from other species to investigate the evolutionary relationship of the Notch protein family. As shown in Fig. 9B, the Notch from insects forms a separate group which includes L. migratoria and O. furnacalis Notch, and the Notch from vertebrates is clustered into another group. Notch from Ciona, sea urchin, and amphioxus is grouped with vertebrate Notch, however, with a low bootstrap value of 82. Caenorhabditis elegans Notch is out of any group, showing a great differentiation from the other taxa. We analyzed the expression profiles of Notch using qRT-PCR methods. As shown in Fig. 7, the expression profile of Notch was similar to that of Hh in two insects. It kept unchanged from the first instar through the fifth instar stage in L. migratoria, while it increased significantly in O. furnacalis fifth instar larvae (Fig. 7). Genes involved in the wg signaling pathway. Wg (the vertebrate homolog of which is Wnt) is transcriptionally activated by Notch signaling, and Wg signaling pathway plays critical roles in axis patterning, cell fate specification, cell proliferation, and cell migration etc [74]. During the development of Drosophila wing, canonical Wg signaling pathway specifies pattern formation along the dorsal/ventral (D/V) axis while Dpp signaling pathway play this role along the A/P axis [75]. The ligand Wg is the founding member of the Wnt family, and is a secreted lipid-modified signaling glycoprotein that has 350-400 amino acids in length [76]. Wg is expressed at the D/V boundary and forms a stable and long-range gradient by symmetrically diffusing at both sides of the boundary. When Wg binds to a receptor complex consisting of the seven-transmembrane protein Frizzled (Fz) and the single-pass transmembrane protein Arrow (Arr, homologous to murine and human low density lipoprotein (LDL) receptor-related protein 5 or 6 (LRP5/6)), the downstream cytoplasmic protein disheveled (Dsh in Drosophila and Dvl in vertebrates) is activated [77]. Dsh in turn inhibits glycogen synthase kinase (GSK)-3b in the b-catenin destruction complex, which mainly consists of Axin, GSK-3b, adenomatous polyposis coli (APC) and b-catenin (armadillo (arm) in Drosophila). Consequently, b-catenin accumulates in the cytoplasm, and stabilized b-catenin then translocates into the nucleus and acts together with the transcription factor Pangolin to regulate the transcription of Wg target genes (Fig. S3D) [77]. In this study, we identified 10 and 11 unigenes for the known components in the Wg signaling pathway from L. migratoria and O. furnacalis transcriptome, respectively (Table 1). No putative Frizzled 2 (fz2) ortholog was identified in L. migratoria transcriptome. A possible reason is that fz2 gene is missing in L. migratoria because of the evolutionary event. The other reason with higher possibility is that the transcript level of fz2 is low and it is not captured in the RNA-seq. We attempted to perform further analysis for identified wg gene, a40139 in L. migratoria and Unigene1219 in O. furnacalis. However, unigene1219 only encodes a 97-amino acid peptide which is too short to have no common sequences with other Wg during the alignment. Therefore, we failed to conduct phylogenetic analysis to reveal the evolutionary relationship of Wg. Additionally, we identified another 12 and 9 unigenes from L. migratoria and O. furnacalis transcriptome, respectively (Table 1), which were potentially involved in the wing development, including homologs to Drosophila apterous, engrailed, homothorax, teashirt, epidermal growth factor receptor, rhomboid, nubbin, pannier, notum, fat, four-jointed, dally-like.

Conclusions
In summary, we sequenced and characterized the transcriptome from the wing discs of L. migratoria nymph. The assembled sequence data comprising 91,907 unique transcripts provides a comprehensive sequence source for future L. migratoria study. We identified a large set of genes relevant to wing development with high significance, especially the genes involved in four signaling pathways -Notch, Hh, Dpp, and Wg signaling pathways, from L. migratoria transcriptome and another O. furnacalis transcriptome obtained previously. The explored wing development-related genes constitute an integrated picture of the development network, which provides the valuable clues for a better understanding of the wing development in L. migratoria and O. furnacalis. These development repertoire genes appear to be evolutionarily conserved to different extent. Functional analyses are necessary to verify our predictions. Nevertheless, the framework of information presented in this study should help to further   The domain organization was predicted using the SMART program (http://smart.embl.de/). Question mark means the end is incomplete. (B) Phylogenetic analysis of the Notch from insects, hydrobios and vertebrates. The unrooted tree was generated using the distance-based method Neighbor-Joing, with Caenorhabditis elegans protein as outgroup. The branch lengths reflect evolutionary divergence. The circled bootstrap value indicates that L. migratoria a1269 and O. furnacalis Unigene6905 (marked in red) belong to insect Notch group. doi:10.1371/journal.pone.0106770.g009 target genes. (C) Notch receptor is a transmembrane protein, with a large extracellular domain (NECD), a transmembrane domain and an intracellular domain (NICD). When its ligands Dl or Ser bind to the Notch, two proteolytic cleavage events are activated: the first cleavage is catalyzed by ADAM-family TACE metalloproteases; the second cleavage is mediated by c-secretase which is an enzyme complex containing presenilin, nicastrin, PEN2 and APH1. The second cleavage liberates the NICD, which then migrates into the nucleus and cooperates with the DNA-binding protein CSL and its co-activator Mam to promote the transcription of downstream target genes. (D) Wg binds to the receptors Fz and Arr. The Wg-Arr-Fz complex binds to and activates the protein Dsh which then disrups the complex of APC/Axin/Sgg and results in an increase of free cytosolic Arm. Arm then translocates to the nucleus where it binds to transcription factor Pangolin (dTCF) and activates the expression of target genes. (EPS)