==== Front
J Insect SciJ. Insect ScijisJournal of Insect Science1536-2442Oxford University Press US 10.1093/jisesa/iey074iey074Research ArticleDe Novo Assembly and Analysis of the White-Backed Planthopper (Sogatella furcifera) Transcriptome http://orcid.org/0000-0002-2839-1774Liang An-Wen 1Zhang Han 2Lin Jia 1Wang Fang-Hai 11 State Key Laboratory for Biocontrol and Institute of Entomology, Sun Yat-sen University, Guangzhou, Guangdong, China2 State Key Laboratory for Biocontrol, Guangdong Key Laboratory of Pharmaceutical Functional Genes, Sun Yat-sen University, Guangzhou, Guangdong, ChinaCorresponding author, e-mail: lsswfh@mail.sysu.edu.cn7 2018 01 8 2018 01 8 2018 18 4 1129 3 2018 © The Author(s) 2018. Published by Oxford University Press on behalf of Entomological Society of America.2018This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.comAbstract
The white-backed planthopper Sogatella furcifera (Horváth) has become an important pest on rice in China and Southeast Asian countries. White-backed planthopper wing bud length is in relation to adult wing length, but little is known about the development and differentiation of wing buds at the molecular level. Using Illumina HiSeq high-throughput sequencing technology, we sequenced four cDNA libraries, two biological replicates of long-winged female fifth-instar nymphs (LW), and two of short-winged nymphs (SW). In total, 62,154 unigenes with an average length of 984 bp and N50 length of 1,878 bp were obtained by de novo transcriptome assembly. A total of 18,416 open reading frames (ORFs) were predicted based on the unigenes. Ninety-three percentage of these ORFs could be annotated by searching for homology in six protein databases. A total of 184 differentially expressed genes (DEGs) with 129 upregulated and 55 downregulated were found in SW compared to LW. Gene Ontology and euKaryotic Orthologous Group classification provided comprehensive information about the function of each gene. Kyoto Encyclopedia of Genes and Genomes enrichment analysis revealed five enriched pathways including three metabolic pathways. In addition, we found that some DEGs were relevant to muscle movement and cuticle and likely involved in development and differentiation of wing buds. This study provided transcriptome resource of female fifth-instar nymphs of white-backed planthopper including long-winged and short-winged nymphs, and different molecular features between them lay the foundation for adult wing morph prediction, promoting further studies on planthopper population management.

illumina sequencingtranscriptomewing budsdifferentiationNational Natural Science Foundation of China10.13039/50110000180931171844Natural Science Foundation of Guangdong Province10.13039/5011000034532015A030313114
==== Body
Sogatella furcifera (white-backed planthopper) has become an important pest on rice in China and Southeast Asian countries (Liu et al. 2000) due to its potential destruction by sucking rice phloem sap and transmitting plant viruses such as ragged stunt virus, grassy stunt virus, and Southern rice black-streaked dwarf virus (Xu et al. 2015a). Similar to another insect in the order Hemiptera, namely Nilaparvata lugens (brown planthopper), the adult white-backed planthopper usually show dimorphism in relation to wing development, occurring in both females and males, though male adults are almost always long-winged (Kisimoto 1956, 1965). The long-winged adults are good at migrating and tracking the growth condition of host plants to avoid the harsh environment. However, short-winged morphs lack of fully developed wings but have a strong ability for reproduction (Liang et al. 2016). Numerous studies revealed that genetic background (Mahmud 1980, Iwanaga et al. 1985, Mori and Nakasuji 1990, Tojo 1991, Yu et al. 1997, Peng et al. 2012) and environmental factors such as temperature (Zhang 1983) often play important roles in planthopper wing dimorphism, but the exact mechanism about wing morph determination has not been illuminated. More recently, Xu et al. (2015b) identified two insulin receptors (insulin receptor 1 [InR1] and 2 [InR2]) to be of great significant to the regulation of wing morph switching in planthoppers. They proposed a model for the molecular regulation of wing polyphenism in planthoppers. The long-winged morph is the basic wing pattern of the winged insects. There are two InRs in the wing buds of rice planthoppers. When the insulin/insulin-like growth factor–like peptide Ilp3 which is secreted by the brain binds to InR1, the insulin signaling pathway is activated leading to long-winged adults. However, when InR2 combined with InR1 to form a heterodimeric receptor, the NlInR1–NlPI(3)K–NlAkt signaling pathway is blocked, thus producing short-winged adults. This model provided a molecular basis for the regulation of wing dimorphism in rice planthoppers, including white-backed planthopper.

In recent years, high-throughput RNA sequencing has greatly facilitated transcriptomics research on insects, especially for insects without reference genome sequences. From 2008 to 2013, insect transcriptomes of 68 species belonging to seven orders have been sequenced by using next-generation sequencing (Zhang and Yuan 2013). Early in this year, the first assembled and annotated whole genome sequence and transcriptome of white-backed planthopper have been reported in public (Wang et al. 2017).

White-backed planthopper wing bud length is in relation to adult wing length (Cook and Perfect 1982, Yamada 1990, Qi et al. 1994), for example, the fifth-instar nymph with wing bud longer than 1.10 mm emerged as long-winged adult, with wing bud shorter than 0.95 mm emerge as short-winged adult (Qi et al. 1994), but little is known about the development and differentiation of wing buds at the molecular level. In this study, we performed de novo transcriptome sequencing on long-winged and short-winged female fifth-instar nymphs of white-backed planthopper by Illumina HiSeq high-throughput sequencing technology. Then raw data were assembled and annotated successively to obtain complete transcriptome information. To identify genes that may be involved in the regulation of white-backed planthopper wing bud differentiation, we performed differentially expressed gene (DEG) analysis between long-winged and short-winged nymphs. These results will be helpful to understand the molecular underpinnings of wing bud differentiation.

As far as we know, it is the first time that the transcriptome of female fifth-instar nymphs of white-backed planthopper including long-winged and short-winged nymphs has been reported. The results in this study will be good resources to learn more about white-backed planthopper physiology and deepen understanding of wing bud development and differentiation in planthoppers.

Materials and Methods
White-backed planthopper Sample Collection and Processing
The white-backed planthopper strain was collected from the rice field located in the South China Agricultural University, Guangzhou, China. Successive generations were reared on rice seedlings under the conditions (28 ± 2°C, 16:8 h light:dark cycle).

One newly emerged long-winged male and one newly emerged short-winged female were mated and cultured in a cage. When the second generation appeared, we selected the 2-d-old fifth-instar female nymphs with the wing buds extended to the fourth uromere more than 0.05 mm as the long-winged nymphs group samples (designed as ‘LW’), but the nymphs’ wing buds were 0.05 mm or more shorter than the fourth uromere as the short-winged nymphs samples (designed as ‘SW’). Each group samples included six individuals, and two biological replicates were set up. All samples were preserved at −80°C prior to RNA extraction.

RNA Isolation and cDNA Library Construction
We extracted total RNA using Trizol reagent according to manufacturer’s instructions (Invitrogen, Carlsbad, CA). Degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Fluorometer (Life Technologies, CA). The integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA). Equal volumes of RNA from each of the two replications of LW and SW were pooled. The mRNA was enriched with magnetic beads with Oligo dT, then fragmented into short fragments. The first-strand cDNA was synthesized by random hexamer primers with the mRNA fragments as templates. With the addition of buffer, dNTPs, RNaseH, and DNA polymerase I, the second-strand cDNA was synthesized and purified with QiaQuick PCR kit and washed off with EB buffer. The purified double-stranded cDNA was subjected to end repair, adding A tail and adaptors, then the products were purified by agarose gel electrophoresis. To obtain the cDNA library, PCR amplification was conducted with the cDNA fragments as templates. Sequencing was done on an Illumina HiSeq 2000 platform.

Data Processing and Analysis
Transcriptome Assembly and Annotation
To ensure the quality of the clean reads for subsequent analysis, reads with adaptor contamination, low quality, and ambiguous base ‘N’ larger than 5% were removed prior to assembly. Q30 and GC content of the clean data were calculated. De novo assembly of short reads without a reference genome was accomplished using Trinity (Grabherr et al. 2011). Open reading frames (ORFs) of all unigenes were predicted by TransDecoder. The ORF and corresponding protein sequence were extracted from the assembled sequences according to the aligned region of the best hit. All the downstream functional analyses were based on the ORFs. Then all ORFs were used as queries to align against sequences in the databases of NCBI non-redundant protein database (Nr), Swiss-Prot, euKaryotic Orthologous Group (KOG), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Ontology (GO) using the BLAST algorithm with a cutoff E value of <10−5. Genes with GO annotations were classified into GO terms, and Web Gene Ontology Annotation Plot (WEGO) (Ye et al. 2006) was applied for plotting GO annotation results.

Identification of DEGs
Bowtie2 (version 2.2.3; Langmead et al. 2009) was employed to align the clean reads used for assembly to the transcript sequences. Then gene expression levels were estimated by RSEM (Li and Dewey 2011) for each sample with a rigorous algorithm called Fragments Per Kilobase of transcript per Million mapped reads (FPKM) which could reduce the effect resulting from differences in gene length and sequencing depth. Differential expression analysis of SW and LW was performed using edgeR (https://bioconductor.org/packages/release/bioc/html/edgeR.html). Unigenes with a false discovery rate (FDR) <0.05 and an absolute value of log2 fold change (FC)≧1 were defined as statistically significant DEGs.

GO and KEGG Enrichment Analysis of DEGs
GO enrichment analysis of DEGs was conducted by using two R packages (GSEABase and GOstats). With regard to KEGG pathway enrichment analysis, R packages—GSEABase and KEGG.db—were adopted. GO term/pathway with FDR <0.05 and Q value <0.05 was selected and defined as significantly enriched GO term/pathway.

Results
Sequencing and De Novo Assembly of the white-backed planthopper Transcriptomes
In order to improve the reliability of the data, two biological replicates for LW and SW were designed for transcriptome sequencing. A total of 163,781,508 high-quality reads were generated for four libraries. Then Q30 was employed to evaluate the base quality. In this study, the average Q30 value of these libraries exceeded 95%, indicating the sequencing data were effective and reliable. The GC content of the nucleotides ranged from 34 to 37%. Then high-quality reads were assembled into 98,097 transcripts and obtain 62,154 unigenes by reducing redundancy. As shown in Table 1, the N50 and mean length of the unigenes were 1,878 and 984 bp, respectively. Unigenes that were 200–400 bp in size accounted for 41.50% of the total number, whereas those longer than 4,000 bp only comprised 3.37%. A total of 18,416 ORFs were predicted based on the unigenes.

Table 1. Summary of the transcriptome assembly

Length range (bp)	Transcripts	Unigenes	
200–400	43,470 (44.32%)	25,794 (41.50%)	
401–600	15,942 (16.25%)	9,703 (15.61%)	
601–1,000	13,361 (13.62%)	8,980 (14.45%)	
1,001–2,000	13,552 (13.82%)	9,294 (14.95%)	
2,001–4,000	8,982 (9.16%)	6,289 (10.12%)	
>4,000	2,780 (2.83%)	2,094 (3.37%)	
Total number	98,097	62,154	
Total length	41,870,291	55,567,878	
N50 length	1,736	1,878	
Mean length	915	984	
Functional Annotation of Transcriptome Sequences
The generated ORFs were annotated by searching and comparing in several databases. As a result, 17,145 (93.1%), 12,808 (69.5%), 12,415 (67.4%), 8,796 (47.8%), 11,436 (62.1%), 9,391 (51.0%) had hits in Nr, Swiss-Prot, GO, Clusters of Orthologous Groups (COG), KOG, and KEGG database, respectively (Table 2). About 60.8% of the ORFs were annotated with sequences from N. lugens; however, there was only 112 ORFs (0.7%) matched to S. furcifera (Fig. 1A). About 44% of the ORFs showed high similarity (80–100%) with sequences in Nr (non-redundant) database (Fig. 1B). E value was used to check the possibility of random matches between two sequences. According to the pie chart in Fig. 1C and E, value ranging from 0 to 1e-100 accounted for more than half.

Table 2. Annotation results of the transcriptome

Database	Annotated number	%	
Nr	17,145	93.1	
Swiss-Prot	12,808	69.5	
GO	184	1.0	
COG	8,796	47.8	
KOG	11,436	62.1	
KEGG	9,391	51.0	
Fig. 1. Characteristics of BLAST searches in Nr database. (A) The species distribution of the Nr annotation results. (B) The similarity distribution of the Nr annotation results. (C) The E value distribution of the Nr annotation results.

To obtain the overview of the functions of annotated ORFs, GO terms of level 2 were assigned to them for classification. It was common that the same ORF got more than one GO terms. Therefore, 43,378 ORFs were categorized into cellular component, 51,526 into biological process, and 15,810 into molecular function. Among these functional terms, cell, cell part, cellular possess, organelle, binding, and metabolic process were six dominant groups (Fig. 2). Similarly, 2,803 ORFs were mapped to two or more KEGG pathways, respectively, whereas 3,733 ORFs were unrelated to any pathway. In addition, 359 KEGG pathways were identified, whose annotated ORF number of more than 100 comprised 13.6%.

Fig. 2. Histogram of GO classification. Three main categories, namely biological process, cellular component, and molecular function, were assigned to 12,415 annotated ORFs.

DEG Identification and Functional Enrichment Analysis
To identify genes involved in wing bud differentiation of white-backed planthopper, we performed a DEG analysis. In total, 184 DEGs were identified, and 129 of them showed a higher expression level in SW compared to LW. FDR <0.05 and an absolute value of log2 FC ≧1 were used as the threshold for selecting the unigenes with significant expression differences. Then we made a volcano plot to intuitively observe the extent of the gene expression difference between two groups and corresponding statistical significance (Fig. 3).

Fig. 3. Volcano plot of DEGs between SW and LW (SW vs LW). The x-axis indicates log2 FC between the two groups and the y-axis indicates the −log10 FDR of gene expression variation. The upregulated genes are shown as red dots, the downregulated genes are shown as green dots and the normal genes are shown as blue dots.

According to GO enrichment results, 19 subfamilies of three major categories were considered to be DEG-enriched terms, which consisted of 12 cellular component terms, 4 molecular function terms, and 3 biological process terms. Nineteen enriched GO terms with Q values and the number of genes in the relevant GO term were displayed in Fig. 4. Among the subcategories, structural molecule activity (26 DEGs) and structural constituent of cuticle (22 DEGs) under molecular function were the two largest groups with the smallest Q values, which suggested that these two GO functions had essential roles in wing buds differentiation of white-backed planthopper. In contrast, virion part, viral capsid, and virion under cellular component had the same and smallest number of DEGs (three DEGs). RNA replication, viral RNA genome replication, and viral genome replication together made up the smallest category, namely, biological process category.

Fig. 4. GO enrichment of the DEGs between SW and LW (SW vs LW). The y-axis indicates the name of GO term and the x-axis indicates rich ratio corresponding to the GO term. The size of the Q value is indicated by the color of the dot. Blue color represents large Q value; red color represents small Q value. The number of DEGs contained in each GO term is represented by the size of the dot.

The KOG database for eukaryotes divided homologous genes from different species into different orthologous clusters based on gene orthologous relationships and with a combination of evolutionary relationship. Consequently, 131 DEGs had hits in KOG database, but only 69 had KOG ID and clustered into 22 function classes. As shown in Fig. 5, R (15.4%), I (11.5%), and K (10.3%) were three major classes.

Fig. 5. Histogram of the KOG classification of the DEGs between SW and LW (SW vs LW). The capital letters in x-axis indicate the KOG categories as listed on the right and the y-axis indicates the number of DEGs in each category.

For KEGG pathway enrichment analysis, only five pathways were significantly enriched and nine genes including three upregulated genes and six downregulated genes were involved in these pathways. The results are displayed in Table 3. In particular, there was one gene enriched in both quorum sensing and peroxisome proliferator-activated receptor (PPAR) signaling pathway. The DEGs involved in quorum sensing and caffeine metabolism were all downregulated in the comparison of SW versus LW.

Table 3. Significantly enriched pathways of the DEGs between SW and LW (SW vs LW)

Pathway	
Q value	DEGs	All genes	
Quorum sensing	0.0157	2 (6.90%)	16 (0.28%)	
Herpes simplex infection	0.0379	3 (10.34%)	118 (2.05%)	
PPAR signaling pathway	0.0379	2 (6.90%)	52 (0.90%)	
Glycine, serine, and threonine metabolism	0.0410	2 (6.90%)	57 (0.99%)	
Caffeine metabolism	0.0451	1 (3.45%)	7 (0.12%)	
Discussion
White-backed planthopper is a small insect in Hemiptera but it has been a huge threat to the rice crops all over the world. To adapt variation in environment, white-backed planthoppers evolve into two wing morphs, and it is well recognized that long-winged adults can fly over a long distance, whereas short-winged adults have robust fecundity. Hence, white-backed planthopper can easily keep away from hash environment and propagate with extreme rapidity to maintain the prosperity of the population. Based on these characteristics, the increasing short-winged adults indicate the expanding of white-backed planthopper population, while the emergence or the increasing of long-winged adults indicates the emigration or immigration phenomenon. Therefore, accurately understanding the wing morph dynamics of white-backed planthopper is beneficial to predict the quantity changes of white-backed planthopper in time and avoid the outbreak of pests. Cook and Perfect (1982) have demonstrated stable relations between the wing bud length of fifth-instar nymph and wing morph of brown planthopper reared on different rice conditions (Qi et al. 1994). Yamada (1990) pointed out that 0.94 mm is the critical value of wing bud length at which fifth-instar nymphs of brown planthopper developed into long-winged female adults. Here, we used RNA-seq to identity genes related to wing bud development and differentiation in white-backed planthopper and then found out the characteristics of these genes by bioinformatics analysis.

In this study, a de novo transcriptome was assembled with sequences from long-winged and short-winged female fifth-instar nymphs. All clean reads were assembled into 98,097 transcripts and 62,154 unigenes of white-backed planthopper. Based on the unigenes, 18,416 ORFs were predicted, of which 17,145 were annotated in Nr database. A comparison with homologous sequences showed that brown planthopper shared the highest similarity with white-backed planthopper in the BLASTx annotation, whereas white-backed planthopper itself had a best match percentage of less than 1%. The reason for this is that brown planthopper has served as a model organism in planthopper research and more sequence resources are available than those of white-backed planthopper in Nr database. According to the KOG classification and KEGG pathway enrichment analysis, it was noted that some DEG sets were involved in transport and metabolism of amino acid and lipid, which meant energy was an important factor for wing bud differentiation. Specially, PPARs are thought to serve as a transcriptional activator of fatty acid catabolism (Contreras et al. 2013). A previous study (Xue et al. 2010) in the transcriptome of brown planthopper has also revealed that PPAR pathway could provide energy for flying by regulating lipid metabolism in macropterous female adults, which played a key role in migration. Later, Contreras et al. (2013) and Yao et al. (2016) found that PPARs also had important effects on glucose and amino acid metabolism. Furthermore, PPARδ can induce genes for long-chain fatty acid oxidation during fasting and endurance exercise in skeletal muscle (Yao et al. 2016), which might be important for regulating muscle activity of long-winged planthoppers in their flight.

Based on the gene expression profiles and Nr annotation, five genes related to muscle movement showed higher expression levels in SW than in LW (Table 4). Among them, three genes showed the most significant upregulated levels in SW. At the same time, they exclusively expressed in short-winged nymphs, suggesting that these three genes were vital genes in wing bud differentiation.

Table 4. DEGs related to muscle movement based on the Nr annotation analysis

Gene ID	Up/downregulation	FPKM in SW	FPKM in LW	FDR	Gene description	
c53337_c4_g1_i4	Up	15.40	0	4.18E-06	Myosin-VIIA isoform X1	
c51535_c1_g1_i1	Up	166.70	4.19	6.17E-03	Troponin C, isotype gamma-like	
c55508_c0_g2_i1	Up	207.97	12.31	3.99E-02	Myosin-2 essential light chain isoform X1	
c54722_c0_g1_i3	Up	56.73	0	3.83E-05	Pro-resilin-like isoform X2	
c44469_c0_g1_i4	Up	65.79	0	1.11E-02	Pro-resilin-like	
In addition, 25 genes with 24 upregulated and only 1 downregulated, encoding cuticle or cuticular proteins, were enriched in SW. All of these genes exhibited dramatically differential expression between SW and LW and three of which expressed uniquely in SW. GO and KOG classification analysis revealed that none of them obtained KOG classification while 18 genes got the same GO term—structural constituent of cuticle in the molecular function category. Cuticle or cuticular proteins had been found in a wide variety of insects, such as Ostrinia furnacalis (Zhang et al. 2016), Athetis lepigone (Li et al. 2013), and Dendrolimus punctatus (Yang et al. 2016) in the order Lepidoptera. These researches all pointed out that many cuticle or cuticular proteins were identified as DEGs in the pairwise comparisons of each developmental stage, highlighting their importance for the development and metamorphosis in insect life history. Both wing buds and cuticle develop from ectoderm. During molting stage, the growth of wing buds coordinated with the metabolism of cuticle (Du et al. 1998). So we hypothesized that cuticle composition differed between SW and LW. The main component of insect cuticle is cuticular proteins and chitin, they together support and maintain the physical structure of the organism and serve as natural barriers against harmful substance from outside world (Moussian 2010, Noh et al. 2016). Previous study showed that poisoned by flufenoxuron, the synthesis of chitin was greatly inhibited (Wang et al. 1996), causing molting disordered and wing buds could not develop into complete and normal wings. In conclusion, there must be close relationships between cuticle and wing buds. Moreover, there were still many crucial genes participating in development of wing buds, including genes associated with sensory process, DNA methylation, and circadian clock.

Conclusion
We assembled the transcriptome of S. furcifera and performed DEG analysis between long-winged and short-winged female fifth-instar nymphs. Six database annotation results helped us comprehensively understand the function of each gene. We finally identified a large number of candidate genes potentially participating in development and differentiation of wing buds, including genes related to muscle movement and cuticle. These findings will provide critical clues for exploring molecular mechanisms of development and differentiation of wing buds. Nymphs at fifth-instar stage are preparing for eclosion to adults. As stated above, the wing bud length of fifth-instar nymphs has complex relationships with the subsequent wing morph of the adults being produced so more research is needed to be done in this area. We hope from some molecular features of wing buds that adult wing morph can be predicted accurately and easily, providing a valuable tool in the study of delphacid population dynamics, thus facilitating the development of pest control.

Acknowledgments
This research was supported by the National Natural Science Foundation of China (31171844) and Natural Science Foundation of Guangdong Province (2015A030313114).
==== Refs
References Cited

Contreras A. V., N. Torres , and Tovar A. R.  
2013 
PPAR-α as a key nutritional and environmental sensor for metabolic adaptation . Adv. Nutr . 4 : 439 –452 .23858092 

Cook A. G. , and Perfect T. J.  
1982 
Determining the wing-morph of adult Nilaparvata lugens and Sogatella furcifera from morphometric measurements on the fifth-instar nymphs . Entomol. Exp. Appl . 31 : 159 –164 .

Du J. G. , Wang Q. , and Cheng X. N.  
1998 
Histological studies on the wing bud development of the brown planthopper Nilaparvata lugens Stål . J. Nanjing Agric. Univ . 21 : 121 –124 .

Grabherr M. G., B. J. Haas M. Yassour J. Z. Levin D. A. Thompson I. Amit X. Adiconis L. Fan R. Raychowdhury Q. Zeng  , et al2011 
Full-length transcriptome assembly from RNA-Seq data without a reference genome . Nat. Biotechnol . 29 : 644 –652 .21572440 

Iwanaga K., S. Tojo T. Nagata  
1985 
Immigration of the brown planthopper, Nilaparvata lugens, exhibiting various responses to density in relation to wing morphism . Entomol. Exp. Appl . 38: 101 –108 .

Kisimoto R  
1956 
Effect of crowding during the larval period on the determination of the wing-form of an adult plant-hopper . Nature 178 : 641 –642 .

Kisimoto R  
1965 
Studies on the polymorphism and its role playing in the population growth of the brown planthopper, Nilaparvata lugens Stål . Bull. Shikoku Agric. Exp. Sta . 13 : 1 –106 .

Langmead B. , Trapnell C. , Pop M. , and Salzberg S. L.  
2009 
Ultrafast and memory-efficient alignment of short DNA sequences to the human genome . Genome Biol . 10 : R25.1 –25.10 .19261174 

Li B. , and Dewey C. N.  
2011 
RSEM: accurate transcript quantification from RNA Seq data with or without a reference genome . BMC Bioinformatics 12 : 1 –16 .21199577 

Li L. T. , Zhu Y. B. , Ma J. F. , Li Z. Y. , and Dong Z. P.  
2013 
An analysis of the Athetis lepigone transcriptome from four developmental stages . PLoS One 8 : 1 –8 .

Liang Z. Q. , Song S. Y. , Liang S. K. , and Wang F. H.  
2016 
Analysis of differential proteins in two wing-type females of Sogatella furcifera (Hemiptera: Delphacidae) . J. Insect Sci . 16 : 1 –7 .26798137 

Liu G. S. , Hanchuan Y. C. , and Shen L. L.  
2000 
Study on wing-morph differentiation of the brown planthopper . Entomol. Knowledge 37 : 186 –190 .

Mahmud F. S  
1980 
Alary polymorphism in the small brown planthopper Laodelphax striatellus (Homoptera: Delphacidae) . Entomol. Exp. Appl . 28 : 47 –53 .

Mori K. , and Nakasuji F.  
1990 
Genetic analysis of the wing-form determination of the small brown planthopper, Laodelphax striatellus (Hemiptera: Delphacidae) . Res. Popul. Ecol . 32 : 279 –287 .

Moussian B  
2010 
Recent advances in understanding mechanisms of insect cuticle differentiation . Insect Biochem. Mol. Biol . 40 : 363 –375 .20347980 

Noh M. Y., S. Muthukrishnan K. J. Kramer , and Arakane Y.  
2016 
Cuticle formation and pigmentation in beetles . Curr. Opin. Insect Sci . 17 : 1 –9 .27720067 

Peng J. , Zhang C. , An Z. F. , Yu J. L. , and Liu X. D.  
2012 
Genetic analysis of wing-form determination in three species of rice planthoppers (Hemiptera: Delphacidae) . Acta Entomol. Sin . 55 : 971 –980 .

Qi L. Z. , Ding J. H. , Ding Z. Z. , Tan H. X. , and Guo Y. J.  
1994 
Researches on relationship between wing bud length of 5th-instar nymph and wing form of adult in Nilaparvata lugens (Stål) . Entomol. Knowledge 31 : 193 –195 .

Tojo S  
1991 
Genetic background of insect migration , pp. 21 –27 . InProceedings, the International Seminar on Migration and Dispersal of Agricultural Insects , National Institute of Agro-Environmental Science , Tsukuba, Japan .

Wang G. Q. , Yan Y. H. , Zhang L. , Zheng Y. L. , and Zhao P.  
1996 
Effect of insect growth regulators on the weight and cuticle of migratory locust . Modern. Agric . 8 : 5 –7 .

Wang L., N. Tang X. Gao Z. Chang L. Zhang G. Zhou D. Guo Z. Zeng W. Li I. A. Akinyemi  , et al2017 
Genome sequence of a rice pest, the white-backed planthopper (Sogatella furcifera) . Gigascience 6 : 1 –9 .

Xu A. L. , Chen S. , He W. A. , Zhang M. Q. , and Qiu Y. F.  
2015a 
Progress on rice resistance to brown planthopper and white-backed planthopper . J. South. Agric . 46 : 1628 –1635 .

Xu H. J., J. Xue B. Lu X. C. Zhang J. C. Zhuo S. F. He X. F. Ma Y. Q. Jiang H. W. Fan J. Y. Xu  , et al2015b 
Two insulin receptors determine alternative wing morphs in planthoppers . Nature 519 : 464 –467 .25799997 

Xue J., Y. Y. Bao B. L. Li Y. B. Cheng Z. Y. Peng H. Liu H. J. Xu Z. R. Zhu Y. G. Lou J. A. Cheng  , et al2010 
Transcriptome analysis of the brown planthopper Nilaparvata lugens . PLoS One 5 : e14233 .21151909 

Yamada S  
1990 
The relation of adult wing-form in the brown planthopper, Nilaparvata lugens Stål (Homoptera: Delphacidae) to wing-pad length of last instar nymphs . Appl. Ent. Zool . 25 : 439 –446 .

Yao Y., Z. B. Wang , and Tang C. K.  
2016 
[PPARs mediate the regulation of energy metabolism by long-chain fatty acids] . Sheng Li Ke Xue Jin Zhan 47 : 1 –6 .27424398 

Yang C. H., P. C. Yang J. Li F. Yang , and Zhang A. B.  
2016 
Transcriptome characterization of Dendrolimus punctatus and expression profiles at different developmental stages . PLoS One 11 : e0161667 .27560151 

Ye J., L. Fang H. Zheng Y. Zhang J. Chen Z. Zhang J. Wang S. Li R. Li L. Bolund  , et al2006 
WEGO: a web tool for plotting GO annotations . Nucleic Acids Res . 34 : W293 –W297 .16845012 

Yu X. P. , X Z. . Lu G. R. Wu , and Tao L. Y.  
1997 
Studies on the immigration and the wing-form development of brown planthopper Nilaparvata lugens (Stål) . Acta Entomol. Sin . 40 : 128 –134 .

Zhang Z. Q  
1983 
A study on the development of wing dimorphism in the rice brown planthopper, Nilaparvata lugens (Stål) . Acta Entomol. Sin . 26 : 260 –267 .

Zhang Q. L. , and Yuan M. L.  
2013 
Progress in insect transcriptomics based on the next-generation sequencing technique . Acta Entomol. Sin . 56 : 1489 –1508 .

Zhang T. , He K. , and Wang Z.  
2016 
Transcriptome comparison analysis of Ostrinia furnacalis in four developmental stages . Sci. Rep . 6 : 1 –9 .28442746