De Novo Assembly and Analysis of the White-Backed Planthopper (Sogatella furcifera) Transcriptome.

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.

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.