Transcriptomic Analysis Suggests Genes Expressed Stage-Independently and Stage–Dependently Modulating the Wing Dimorphism of the Brown Planthopper

Wing dimorphism is considered as an adaptive trait of insects. Brown planthoppers (BPHs) Nilaparvata lugens, a serious pest of rice, are either macropterous or brachypterous. Genetic and environmental factors are both likely to control wing morph determination in BPHs, but the hereditary law and genes network are still unknown. Here, we investigated changes in gene expression levels between macropterous and brachypterous BPHs by creating artificially bred morphotype lines. The nearly pure-bred strains of macropterous and brachypterous BPHs were established, and their transcriptomes and gene expression levels were compared. Over ten-thousand differentially expressed genes (DEGs) between macropterous and brachypterous strains were found in the egg, nymph, and adult stages, and the three stages shared 6523 DEGs. The regulation of actin cytoskeleton, focal adhesion, tight junction, and adherens junction pathways were consistently enriched with DEGs across the three stages, whereas insulin signaling pathway, metabolic pathways, vascular smooth muscle contraction, platelet activation, oxytocin signaling pathway, sugar metabolism, and glycolysis/gluconeogenesis were significantly enriched by DEGs in a specific stage. Gene expression trend profiles across three stages were different between the two strains. Eggs, nymphs, and adults from the macropterous strain were distinguishable from the brachypterous based on gene expression levels, and genes that were related to wing morphs were differentially expressed between wing strains or strain × stage. A proposed mode based on genes and environments to modulate the wing dimorphism of BPHs was provided.


Introduction
Polyphenism is a life history strategy for organisms to deal with heterogeneous environments. Two or more distinct morphs can arise from a genotype as a result of differing environmental conditions [1]. Wing dimorphism provides an opportunity for insects to choose migration or settlement. An aphid clone can produce both the wingless and winged offspring under different conditions, such as population density [2,3], and a pair of rice planthopper parents can generate two kinds of progeny with long or short wings [4,5]. The waterstrider Limnoporus canaliculatus can produce both the long-winged and wingless morphs that are determined by genetic component and photoperiod [6]. It seems that the wing dimorphism in insects is a phenotypic plasticity of a morphological trait, but it has been confirmed in multiple species that the wing morph is solely determined by genetic mechanisms, or solely by environmental mechanisms, or through a combination of both [7][8][9]. In Myzus persicae, mitochondrial adenine nucleotide translocase (ANT), the chemoreception and takeout-like (TOL) genes showed higher expression levels in the winged morphs than in the wingless morphs [10]. The genetic and molecular mechanisms underlying wing dimorphism in insects are attracting more attention in the last decade, due to the development of molecular techniques [7,10,11]. However, the genetic law and genes network determining wing polymorphism are still unclear.
The rice brown planthopper (BPH), Nilaparvata lugens, which is a devastating pest of rice, displays an obvious wing dimorphism, and it is a better model organism to study on the wing morph determination [7]. In a wild population, the long-and short-winged morphs often coexist [12]. Previous researches have shown that genetic and environmental factors together contribute to the wing morph of rice planthoppers [5,13,14]. The nutrient of host plants, population density, and climate are all factors that are known to affect wing morphs of insects [13][14][15][16][17][18][19][20]. Phenotypic differentiation of wing morphs might be triggered by environment cues occurring in some specific developmental stages of insects [21][22][23][24][25]. Wing induction in aphids can be controlled either by the mother (pre-natal) or by the developing nymph (post-natal) crowding, depending on aphid species [23,26]. The 3rd or 4th-instar nymphs of BPHs are the most sensitive to a decrease of population density, which induces the short-winged morph [22,24]. The determination of wing morphs, long-or short-winged, and winged or wingless in insects might be made in a short window period before adult emergence, and this period is generally named the 'sensitive stage'. These results imply that gene roles in determining wing morphs might be various in different developmental stages of insects. Therefore, it is worth studying changes in the gene expression levels between the long-and short-winged insect strains across all developmental stages to highlight the molecular mechanism underlying wing dimorphism.
Gene expression differences between the winged and wingless, or long-and short-winged morphs have been illustrated in some insects [27,28]. More than one-thousand of differentially expressed genes were identified between two wing morphotypes [27][28][29]. For example, 1663 differentially expressed transcripts were identified from the winged and wingless cotton aphids Aphis gossypii, while using a tag-based digital gene expression approach, and these transcripts were enriched in the metabolic pathways of ribosome, pyruvate metabolism and proteasome, protein synthesis and degradation, lipid metabolism, immunity, RNA transport, and some signaling pathways [28]. In 1734 unique cDNAs (an estimated 10% of coding genes) from the pea aphid Acyrthosiphon pisum, there were 141 and 142 genes with differential transcript accumulation between winged and wingless morphs in the fourth instar nymph and adult, respectively, and the functions of differentially accumulated transcripts mainly enriched in muscles and energy production [27]. Six out of 11 wing development-related genes showed significant differences in expression level among five developmental stages of the pea aphid, and another two genes exhibited a significant development stage × wing morph interactive effect [30]. In the New Zealand stonefly Zelandoperla fenestrata, the fully winged and vestigial-winged morphotypes strongly differentiated in small regions of the genome [31]; nine and one wing development-related gene clusters were significantly upregulated in full-winged and vestigial-winged ecotypes, respectively [32]. In the third-instar nymphs of BPHs, there were 2544 differentially expressed genes between nymphs that were reared on yellow-ripe stage of rice with a high population density (inducing long-winged morphs) and on tillering stage rice with a low population density (inducing short-winged morphs), and the expression levels were dependent on the developmental stage of nymphs [29,33]. The molecular mechanism underlying wing polymorphism in insects is still not clear [34], although a previous study has known that the expression of two insulin receptors (InR1 and InR2) determined wing morphs of rice planthoppers [7].
Interactions between genetics and environment conditions manipulate the wing morph of rice planthoppers [4,18,[35][36][37]. The direct selection for wing morphs in BPHs under a high-density condition was successful in obtaining the long-winged and short-winged pure-bred lines [38]. Therefore, using the pure-bred lines to study on variations in transcriptome and gene expression levels between the two wing morphotypes would be much better than using the wild populations to evaluate the genes network in determining the wing morphs of insects. In this study, we hypothesized that genetic pathways that enriched differentially expressed genes globally or stage-dependently across egg, nymph and adult stages determined the wing dimorphism of BPHs. Therefore, the long-and short-winged pure-bred or nearly pure-bred strains of BPHs were set up while using the successively directional selection, and then transcriptomes of their eggs, 3rd-instar nymphs and adults were sequenced. The differentially expressed genes and their enrichment pathways between the two strains and among the three growth stages were analyzed. The results will provide more evidence for the genetic and molecular mechanisms of wing dimorphism in insects.

Insects
The brown planthoppers (BPHs), Nilaparvata lugens, were collected from rice field in Nanjing, China and reared in a chamber under a 14 h light:10 h dark photoperiod at 25 ± 1 • C and a relative humidity of 75% ± 10% using rice seedlings (rice variety: Wuyunjing-7) [5,14]. The newly born long-winged and short-winged adults were used in the selection experiments for wing morphs.

Selection of Macropterous and Brachypterous Pure-Bred Strains
The long-winged (macropterous) and short-winged (brachypterous) strains were established through sibling inbreeding for successive generations under a constant condition [5]. In each selection generation, 10 pairs of adults and 30 first-instar nymphs from each pair of adults were examined. During a selection process, a pair of newly emerged unmated female and male adults with the same wing form were paired (M♀× M♂or B♀× B♂), copulated and produced progenies in a plastic cup (diameter 80 mm, height 100 mm) with 10 rice seedlings. Ten independent pairs were set up for each selected lineage. Within 24 h after the eggs hatching, 10 first-instar nymphs were collected and reared in a plastic cup with 10 rice seedlings, and three replications of nymphs were performed for each pair of adults. In each generation, the wing form of all emerged adults was examined and the short-winged morph rate was calculated. The offspring adults with an identical wing form as their parents were selected for sibling mating to produce the next generation. After 38 and 36 successive generations of selection in the long-winged and short-winged lineages, respectively, the percentages of the short-winged adults were 13.98% ± 3.33% and 95.67% ± 2.01% ( Figure 1). These long-winged and short-winged lineages were considered as the nearly pure-bred macropterous (M) and brachypterous (B) strains, and their eggs (named ME and BE), 3rd-instar nymphs (M3rd and B3rd), and adults (MA and BA) were used for the transcriptome sequencing and qPCR experiments. Additionally, the selection of wing morphs remained.

RNA Extraction and cDNA Library Construction for Transcriptome Sequencing
The samples used for transcriptome sequencing were collected from the 39th generation of macropterous strain and 37th generation of brachypterous strain. The eggs, third-instar nymphs, and adults from the macropterous and brachypterous strains were examined. Adults from the 38th generation-selected macropterous strain and the 36th generation-selected brachypterous strain were chosen to produce eggs on rice seedlings. The seedlings were renewed at a 24 h interval and dissected seedling tissue to collect eggs five days later, and these fertilized eggs with red eyepoints were dissected out from rice leaf sheaths while using a pair of ophthalmic forceps under a dissection microscope. 900 eggs were collected for extracting RNA. 40 nymphs were collected when the nymphs grew up to the 3rd instar. When adults emerged, the 24-36 h-old adults (15 females and 15 males) were collected and pooled as an adult sample for transcriptome sequencing. All of the samples were frozen in liquid nitrogen and then stored at −80 • C.
Total RNA was isolated while using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and treated with DNase I (Takara, Dalian, China), according to the manufacturer protocol. The purity and integrity of total RNA were determined by spectrophotometer NanoDrop ND-8000 (Thermo Fisher, Waltham, MA, USA) and Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). The samples with RNA integrity number (RIN) > 8.0 were used for sequencing. More than 20 µg RNA for each sample was obtained to construct the cDNA library.

cDNA Library Construction and Sequencing
Six cDNA libraries were separately constructed according to the TruSeq RNA Sample Prep Kit v2 (Illumina, San Diego, CA, USA). 200 ng total RNA sample was purified by oligo-dT magnetic beads, and then poly (A)-containing mRNA were fragmented into small pieces with Elute, Prime, Fragment Mix. The first-strand cDNA was synthesized while using First Strand Mater Mix and Super Script II (Invitrogen, Carlsbad, CA, USA) reverse transcription (The reaction condition: 25 • C for 10 min.; 42 • C for 50 min.; 70 • C for 15 min.), and then added the Second Strand Master Mix to generate the second-strand cDNA (16 • C for 1 h). The purified fragmented cDNA was combined with the End Repair Mix and then incubated at 30 • C for 30 min. The end-repaired DNA was purified with Ampure XP Beads (Agencourt, Danvers, MA, USA) and then added the A-Tailing Mix and mixed well by pipetting and incubated at 37 • C for 30 min. The cDNA fragments with Poly (A) addition was combined with the Adenylate 3 Ends DNA, RNA Index Adapter and Ligation Mix, and then mixed well by pipetting to connect the sequencing adaptor. The ligate reaction was incubated at 30 • C for 10 min. The end-repaired DNA was purified with Ampure XP Beads, and then several rounds of PCR amplification with the PCR Primer Cocktail and PCR Master Mix were performed to enrich the cDNA fragments. The PCR products were purified with Ampure XP Beads and created the final cDNA library. The qualified library was amplified on cBot to generate the cluster on the flowcell (TruSeq PE Cluster Kit V3-cBot-HS, Illumina), and the amplified flowcell was sequenced pair end on the HiSeq2000 System (TruSeq SBS KIT-HS V3, Illumina) at Beijing Genomics Institute (Shenzhen, China). The read length was 90 bp.

De novo Assembly and Annotation
Based on the eggs, third-instar nymphs and adults from both the macropterous and brachypterous strains, 2.8 × 10 10 bases in total were generated in all six transcriptomes. The clean reads were obtained via removing reads with adaptors and unknown nucleotides larger than 5%, and low-quality reads (the rate of bases with quality value ≤ 10 was more than 20%). The proportions of unknown nucleotides in clean reads were 0.01% for all six samples. Transcriptome de novo assembly was carried out while using these short clean reads on Trinity program (version v2.0.6) and attained the unigenes. Subsequently, the unigenes were taken into further process of sequence splicing, redundancy removing and clustering with TGICL tools version 2.1 [39,40]. Unigene sequences were firstly aligned to protein databases NR, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Clusters of Orthologous Groups (COG) (e-value < 0.00001) by Blastx (version v2.2.23), and nucleotide database NT (e-value < 0.00001) by Blastn (version v2.2.23), retrieving proteins with the highest sequence similarity with the given unigenes, along with their protein functional annotations. The sequence direction and the coding regions of unigenes were determined according to the best aligning results. The ESTScan program predicted unigenes unaligned to none of the above databases [41]. The unigenes were mapped to the COG database and predicted the possible functions. With NR annotation, Blast2GO program (version v2.5.0) was used to get Gene Ontology (GO) annotation of unigenes. WEGO software [42,43] was used to undertake GO functional classification for all unigenes and understand the distribution of gene functions of the species from the macro level after obtaining GO annotation for every unigene. With the help of the KEGG database, genes biological complex behaviors were further studied, and using KEGG annotation we obtained pathway annotation for unigenes. InterProScan5 performed the InterPro annotation (version v5.11-51.0).

Analysis of Differential Gene Expression Profiles
Bowtie2 remapped clean reads to the unigenes (version v2.2.5). The relative expression levels of all the matched unigenes were normalized by transforming the clean data to fragments per kilobases of transcripts per million mapped fragments (FPKM) by RSEM (version v1.2.12). The log 2 (fold change) of FPKM of an unigene in one transcriptome to another was used to determine the differentially expressed gene (DEG). The DEGs between the macropterous (M) and brachypterous (B) strains in the egg, 3rd-instar nymph, and adult stages were identified based on the Poisson distribution method [44]. False discovery rate (FDR) was calculated to determine the threshold p-value in multiple test. Transcripts with a minimal two-fold change in expression (|log 2 fold change| ≥ 1) and FDR ≤ 0.001 were considered as DEGs between two samples.
GO and KEGG pathway enrichment of DEGs were performed. This enrichment analysis would find the main biological functions of DEGs, and the main biochemical pathways and signal transduction pathways that DEGs took part in. The enriched p-values were calculated according to the hypergeometric test and performed with Bonferroni correction. The GO terms with corrected p-value < 0.05 were defined as significantly enriched GO terms in DEGs. The REVIGO web tool was used to summarize the long lists of GO terms with the default parameters in order to remove redundant GO terms [45]. For the pathway enrichment, the multiple testing correction Q value < 0.05 was used as the threshold.
The principal component analysis (PCA) was performed for the six samples (ME, M3rd, MA, BE, B3rd, and BA) based on the FPKM of each unigene after data standardization while using the z-score method. The first two principal components (PC1 and PC2) interpreted 64.9% of variances. Therefore, we used the PC1 and PC2 to distinguish the six samples (egg, 3rd-instar nymph, and adult from the macropterous and brachypterous strains). The hierarchical clustering of six samples was analyzed based on the expression levels of 6523 DEGs shared with eggs, nymphs, and adults.
Expression pattern analysis (trend analysis) was performed to obtain the transcriptional differences over the developmental stages between the macropterous and brachypterous strains. All DEGs across three developmental stages were assigned to eight expression profiles while using a short time-series expression miner STEM (version v1.3.8). The unigenes belonging to the same expression profiles had similar expression pattern among the growth stages. For each strain, the clustered profiles with p < 0.05 were considered as the significant trend profiles. GO and KEGG enrichment were performed in the trend profiles.

Expression Levels of Genes Related to Wing Development in Macropterous and Brachypterous Strains
Eight differentially expressed transcripts that were related to the wing development of insects [7,30,46] were selected for measuring their expression levels in the macropterous and brachypterous strains during the egg, 3rd instar nymph, and adult stages while using the qPCR method ( Table 1). The total RNA was extracted from eggs, third-instar nymphs, female adults, and male adults from the macropterous and brachypterous strains. A total of 1000 ng total RNA from each sample was used to synthesize the single-strand cDNA while using PrimeScript RT reagent kit with gDNA Eraser (TAKARA, Dalian, China). 18S rRNA gene was used as an internal reference gene. Gene-specific primers were designed with Beacon Designer (version 7.0) and Oligo (version 7.0) software (Table 1), and the primer specificity was verified while using the dissociation curve analysis. The qPCR was performed using SYBR Premix Ex Taq kit (TAKARA, Dalian, China) according to the manufacture's protocol in an ABI 7500 (Applied Biosystems, Carlsbad, CA, USA). Each reaction mixture was 20 µL containing 10 µL SYBR Premix Ex Taq, 200 nM of each forward and reverse primer, 0.4 µL ROX reference DyeII. The cycling parameters were 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s, 60 • C for 34 s. All of the reactions were performed in triplicate, and dissociation curve analysis was performed after each assay to determine the target specificity. The cycle threshold (Ct) was normalized to the 18S rRNA gene (∆Ct). The Ct values of 18S rRNA gene were stable between the macropterous (14.45 ± 0.20) and brachypterous (14.62 ± 0.22) strains. The ∆Ct of the macropterous sample was calibrated against the corresponding brachypterous sample to obtain the ∆∆Ct value, and this value was used to evaluate the fold-change of gene expression level in the macropterous strain related to the brachypterous strain. The ∆Ct value from a sample was calibrated again with one sample from a specific stage of BPHs and obtained another ∆∆Ct, and the 2 −∆∆Ct was used to calculate the relative expression level of a gene in a sample. Three biological replications were performed for all of the experiments. The correlation of the fold changes of eight genes between the macropterous and brachypterous strains measured by the RNA-seq and qPCR was analyzed while using the Pearson method. The maximum of the fold change from a female sample and a male sample measured by qPCR was considered as the fold change of an adult sample because females and males were pooled in an adult sample for the RNA-seq. The effects of the BPH strain and developmental stage on the fold changes were analyzed while using GLM, and the differences in expression levels of each gene between the macropterous and brachypterous strains were analyzed using student t-test. All of the statistics were performed in IBM SPSS Statistics 25.

Selection Response of Wing Morphs
The selection response of BPHs in wing morphs was strong ( Figure 1). The frequency of the short-winged morph increased up to 92-100% and maintained during 3-60 generations of directional selection in the brachypterous female × brachypterous male lineage. The frequency of the long-winged morph could increase up to approximately 95% after 20 generations of directional selection in the macropterous female × macropterous male lineage, nevertheless the frequency still fluctuated during 29-60 generations (Figure 1).

Differentially Expressed Genes between the Macropterous and Brachypterous Strains
Illumina sequencing generated 54.42-57.60 million of 90 bp pair-end raw reads from the six samples (NCBI SRA accession: SRR10008513), and the total mapped reads were 93.03-94.49%. In the results of assembly, 77,765 unigenes were detected ( Table 2). The total length for unigenes was 102,867,302 nt with 1323 nt of average length and 2482 nt of N50. The total annotation unigenes were 34,589. The wing dimorphism resulted in expression differentiation of many genes in BPHs across three growth stages. As shown in Supplementary Materials Tables S1-S3, 13,606, 14,706, and 15,683, DEGs were identified between the macropterous and brachypterous strains in the eggs, 3rd instar nymphs, and adults, respectively, accounting for 17.5, 18.9, and 20.2 percent of total unigenes (Figure 2).  The total number of DEGs increased as the development of brown planthoppers from egg to adult. In eggs and third instar nymphs, nearly the same numbers of genes were upregulated vs. downregulated in the macropterous strain as compared to the brachypterous strain ( Figure 2). In contrast, substantially more genes were upregulated in macropterous adults as compared to brachypterous adults (Figure 2).

GO Term and KEGG Classification for DEGs between the Macropterous and Brachypterous Strains
The GO term showed that the DEGs between the macropterous and brachypterous strains were distributed across the similar biological process, cellular components, and molecular functions across the three different growth stages of BPHs egg, nymph, and adult ( Figure 3). While using the Revigo web tool, we found that significantly enriched GO terms were mainly involved in cell fate determination, regulation of cell projection organization, site of polarized growth, growth cone, microtubule binding, and phosphate transmembrane transporter activity in eggs, whereas significantly enriched GO terms were the metabolic process, catalytic activity, and oxidoreductase activity, etc. in the 3rd-instar nymphs and adults (Table 3).
Similarly, the numbers of DEGs between the macropterous and brachypterous strains were classified into similar KEGG pathways among the three growth stages of BPHs. Overall, many DEGs were involved in metabolism, environmental signal transduction, bacterial infection, endocrine, and digestive systems ( Figure 4).
The DEGs were significantly enriched in four KEGG pathways: regulation of actin cytoskeleton, focal adhesion, tight junction, and adherens junction across three growth stages of BPHs: egg, nymph, and adult ( Figure 5). Additionally, the insulin signaling pathway and amoebiasis were enriched with DEGs in the egg stage ( Figure 5A), and many DEGs were enriched in metabolic pathways, amoebiasis, vibrio cholera infection, vascular smooth muscle contraction, platelet activation, oxytocin signaling pathway, amino sugar and nucleotide sugar metabolism, and glycolysis/gluconeogenesis in the 3rd-instar nymph stage ( Figure 5B). In the adult stage, the DEGs were enriched in metabolic pathways, vascular smooth muscle contraction, protein processing in endoplasmic reticulum, lysosome, carbon metabolism, and amino sugar and nucleotide sugar metabolism ( Figure 5C).

Molecular Differentiation of Macropterous and Brachypterous Strains
Wing morph selection caused remarkable changes in the gene expression profiles of BPHs. A total of 6523 DEGs were found to consistently differ between the macropterous and brachypterous strains across all three growth stages: egg, 3rd-instar nymph, and adult ( Figure 6A). The six samples from the brachypterous (BE, B3rd, and BA) and macropterous (ME, M3rd, and MA) strains were discriminable based on the first two components that were derived from the FPKM of all 76,742 unigenes (Table S4), and first principle components axis PC1 clearly discriminated based on life stage, whereas the second axis PC2 discriminated based on wing morph ( Figure 6B). The eggs, 3rd-instar nymphs, and adults could be clustered well into the brachypterous or macropterous strain while using the expression levels of 6523 DEGs, suggesting significant differences in gene expression levels between the wing strains ( Figure 6C). The trend analysis for all unigenes among egg, nymph, and adult stages showed that there were three significant trend profiles (profile 0, 1, and 3) out of eight profiles (profile 0-7) in both the macropterous ( Figure 7A) and brachypterous ( Figure 7B) strains. The general significant expression trends for genes were for a decline from the egg to the adult stage ( Figure 7A,B). However, it is important to note that the numbers of genes distributed within each of the three significant trend profiles were different between the macropterous and brachypterous strains. Many genes were distributed to the profile 0 in the brachypterous strain, whereas to profile 3 in the macropterous strain ( Figure 7A,B). Moreover, genes that were distributed to the same trend profile in the brachypterous and macropterous strains were usually involved in different KEGG pathways ( Figure 7C). The expression profiles for DEGs between the macropterous and brachypterous strains did not change across the three growth stages for six pathways: platelet activation, insect hormone biosynthesis, adherens junction, regulation of actin cytoskeleton, insulin signaling pathway, and PI3K-Akt signaling pathway (e.g., they were constitutively different between strains). In contrast, the DEGs that were enriched in six other pathways did show different expression profiles across the three growth stages for the macropterous and brachypterous strains (tight junction, focal adhesion, vascular smooth muscle contraction, glycolysis/gluconeogenesis, oxytocin signaling pathway, and metabolisms). Meanwhile, the genes in yet a third group of pathways were not differentially expressed between the strains while at the same growth stage, but the expression profiles were different when being considered across the three growth stages (e.g., for circadian rhythm-fly, cell adhesion molecules, hippo signaling pathway-fly, Wnt signaling pathway, etc.). For example, genes in the circadian rhythm-fly pathway belonged to profile 3 in the brachypterous strain, but they belonged to profile 0 in the macropterous strain ( Figure 7C).

Expression Levels of Genes in Macropterous and Brachypterous Strains
The qPCR results showed that fold changes in the expression level of eight selected genes measured by RNA-seq method were significantly related to the values measured by qPCR method in three growth stages ( Figure 8A), which suggested the validation of the RNA-seq in this study. The wing-morph strain, growth stage, or the interactions significantly affected the relative expression levels of these six genes chico, InR1, torc1, srf, tyr1, and foxo ( Figure 8B,C,E-G,I). But the relative expression level of hth was not significantly different between the macropterous and brachypterous strains, and among eggs, 3rd instar nymphs, and adults ( Figure 8H). The InR2 was differentially expressed among three growth stages of BPHs, and the wing morph affected the expression of InR2 associated with the growth stage ( Figure 8D). The expression levels of chico and torc1 were higher in the brachypterous strain than that in the macropterous strain during all three growth stages ( Figure 8B,E), and the other genes at least were higher in a specific growth stage (Figure 8C,F,G,I).

Discussion
Differentiation in the wing morph of insects leads to remarkable changes in the gene expression levels [7,30,33,39]. In this study, we used the pure-bred strains of the short-winged and long-winged BPHs to explore the global changes in gene expression, and found that there were substantial and consistent changes between the macropterous and brachypterous strains across three growth stages, egg, nymph, and adult, but there were also substantial growth stage-dependent changes. These changes in gene expression occurred in the genetic pure-bred strains of wing morphs under a constant environmental condition might result from the genetic basis. Here, we found 15,683 DEGs between adults (including females and males) from the macropterous and brachypterous strains that were twenty times more than the 755 DEGs found between the female macropterous and brachypterous adults measured by Xu et al. [33]. We thought that the number of DEGs between the macropterous and brachypterous pure-bred strains might be larger than the wild line, although there also were 2172 genes with significant expression changes between the macropterous female and male adults [33]. The multigenerational selection for wing morphs may induce more genes to express differentially between the macropterous and brachypterous morphs. Wing selection produced lines with a relatively consistent wing form in BPHs [5,47], which showed the importance of genes in determining wing morphs. From this study, many differentially expressed genes were found between the macropterous and brachypterous strains of BPHs, which suggested that gene expression might determine the wing morphs of BPHs. There were 3882 out of 7909 genes expressed differentially in the functional and histolytic thoracic flight muscles of Gryllus firmus, and the fat body of short-winged and long-winged morphs of this insect also exhibited transcriptional differences [48]. These findings can explain well the previous study that the wing morphism of some insects, including rice planthoppers, is under the polygenic control [6,35,49].
Although, in this study, the RNA-seq samples were collected from the 37th and 39th generation of wing selection that were nearly pure-bred lines of the brachypterous and macropterous morphs and the threshold to discriminate DEGs was strict (fold change ≥ 2 and FDR ≤ 0.001), the presence of false-positive DEGs was still not wholly avoid, due to the absence of biological replicates. The high percentages of DEGs in each developmental stage (17-20%) between the brachypterous and macropterous strains may include a few false positive errors. Ideally, the results will be better if it was repeated two additional times with independently inbreed strains, or at least on different generations of the original inbreed strains to account for the variability in the short-winged rate. In this study, we used the FPKM to estimate the expression level of genes, but the Transcripts per Kilobase Million (TPM) might be more reliable than FPKM or RPKM [50]. We reanalyzed the data while using TPM, and DEGs analyzed by TPM were as similar as these by FPKM. A total of 13,606, 14,672, and 15,706 DEGs between the macropterous and brachypterous strains in eggs, 3rd instar nymphs, and adults were found based on TPM, respectively, and it was 13,606, 14,706, and 15,683 based on FPKM (Tables S1-S3). Accordingly, we thought that the methods to estimate the DEGs in this study were suitable and the false-positive error might be low. The method DEseq2, EdgeR, or Limma was recently used for DEGs analysis based on the negative binomial distribution model [32,51], but we performed the analysis with the Poisson distribution method previously used in other studies [52,53], followed the qPCR test for some putative target genes in this study. The old method seemed to be suitable for this RNA-seq data, but the new one might be much better. Therefore, the comparison between different methods of DEG analysis for this RNA-seq data might be considered.
We found in this study that the DEGs between the macropterous and brachypterous strains were mainly distributed to the similar GO terms, such as catalytic activity, binding, metabolic process, and cellular process, and the similar KEGG units, such as global and overview maps, signal transduction, cellular community, endocrine system, and digestive system among the egg, 3rd instar nymph, and adult stages. The results showed that these functional terms and units might be involved in the regulation of wing dimorphism, and these genes related to the cellular and metabolic process, signal transduction, endocrine, and digestion may be the key genetic basis of wing morphs. In aphids, the metabolic and signaling pathways, muscles, and energy production were enriched for DEGs between the winged and wingless morphs [27,28]. Moreover, the significantly enriched GO terms for DEGs between the macropterous and brachypterous BPH strains in eggs were much more than that in nymphs and adults (Table 3). This result showed that genetic differentiations in the molecular level between the macropterous and brachypterous strains were significant in the developmental stage of egg, which suggested the genetic determination of wing morphs in BPHs. Gene profiles of the macropterous and brachypterous strains across developmental stages, to our knowledge, were first identified here, although the genetic basis or inheritance of wing dimorphism in some species of insects, including rice planthoppers, has been found [2,5,35,47,49].
Which are the genes involving in the determinant of wing morphs in insects? This is the key question to address the gene network and molecular mechanism of wing dimorphism. In this study, we found that there might be at least two groups of genes that were differentially expressed between the pure-bred macropterous and brachypterous strains. One group was the constant differentially expressed genes across three growth stages, and here we called 'global genes', which were mainly related to the regulation of actin cytoskeleton, focal adhesion, tight junction, and adherens junction. These genes may determine the morphological structure of wings and muscles via mediating various important cellular processes, such as cell structural support, axonal growth, cell movement, cell adhesion, and so on [54][55][56], and they result in the different wing sizes and flight muscle in the macropterous and brachypterous adults [57,58]. We found that the expression levels of chico and torc1 genes in the brachypterous strain were significantly higher than that in the macropterous strain, regardless of the BPH growth stages. The gene chico acted autonomously in the control of cell size and organ size [59]. RNAi chico in BPHs would result in the short-winged morph [7]. The torc1 gene regulates multiple cellular processes to control cell growth in response to environmental signals [60].
The genetic information in these global genes might represent the maternal determinant for progeny wing morphs.
The other group genes were differentially expressed between the macropterous and brachypterous strains only in a specific growth stage, named here 'stage-dependent genes', which mainly involved in the metabolic pathways. These pathways, the insulin signaling pathway, oxytocin signaling pathway, amino sugar and nucleotide sugar metabolism, carbon metabolism, and glycolysis/gluconeogenesis were significantly enriched with DEGs between the macropterous and brachypterous strains in either eggs or the 3rd instar nymphs. We speculated that the expression of these stage-dependent DEGs might be sensitive to environmental factors, including host quality and crowding, and, consequentially, their expression levels may determine the wing morph. The heritability of wing morphs in many insects was relative low and inconsistent [47,49], which suggested a threshold in the gene expression level might modulate wing forms. A study reported that two insulin receptors (InR1 and InR2) determine the wing morph of the rice planthoppers via regulating the activity of the forkhead transcription factor foxo, and the knockdown of InR1 gene led to the short-winged morph, but the knockdown of InR2 gene led to the long-winged morph [7]. However, the knockdown of foxo directly resulted in the long-winged morph whether the InR1 and InR2 were interfered or not [7]. The srf transcription factor regulates or promotes tissue growth in Drosophila, which is largely based on the overexpression of Pico and Mal reported to increase wing size [61,62]. In this study, we found that the interaction between wing morph strain and developmental stage of BPHs determined the expression levels of InR1, InR2, foxo, and srf genes. Environmental conditions, especially host quality, such as glucose concentration, affected the wing morph of BPHs [19]. These environmental conditions may affect insect's wing morphs via regulating the expression levels of stage-dependent genes. Therefore, environmental changes that only occur in a specific period of insects can result in the wing morph shift between macropterous and brachypterous [21][22][23][24].
The nature of the growth stage-dependent genes versus the globally differentially expressed genes between the macropterous and brachypterous strains suggests that the wing morph of BPHs might be regulated via a two-stage genetic path: the first stage is the inheritance from adults that determines the expression of all 'global genes', and the other is the environmental effect in nymphs that determines the expression of all 'stage-dependent genes' (Figure 9). Among global DEGs, these related to the tight junction, focal adhesion and vascular smooth muscle contraction pathways exhibited different expression trend profiles across the egg, 3rd instar nymph, and adult stages between the macropterous and brachypterous strains, and the others related to the adherens junction and regulation of actin contraction pathways were constantly expressed across the three growth stages of BPH. These global genes may all be intimately involved in the determination of structural foundation of wings, and their expression levels were under genetic control. These growth stage-dependent genes that are related to metabolic pathways, such as insulin signaling pathway, oxytocin signaling pathway, amino sugar and nucleotide sugar metabolism, and carbon metabolism, may finally determine the phenotype of wings via responses to environments in nymphs. We speculate that these genes may be sensitive to environmental conditions only in the heterozygous wing line, but insensitive in the homozygous wing strain, like the global genes. Therefore, the genetic pure-bred line in wing morph can be selected and their offspring's wing form was not altered by different environmental conditions, such as photoperiod [14,38,47,49]. The genetic and environmental factors may modulate the wing morphs of BPHs via cooperation of two types of genes expressed stage-dependently and stage-independently ( Figure 9). In this study, we found that, after 59 generations of wing morph directional selection, there still were an extremely small amount of short-winged and long-winged individuals that were produced by the macropterous and brachypterous strains, respectively ( Figure 1). The result suggests that the pure-bred wing morph are not easy to be fully segregated in BPHs. Mochida (1975) found a BPH strain produced abundant brachypterous adults [63]. Even if the InR1 or InR2 was knocked down, not all of the treated BPHs showed a wholly identical wing morph [7]. In field survey, we found that there often were some BPHs on rice that did not grow up to long-winged morphs in fall to migrate, even though they could not overwinter there. These non-migrant BPHs may be genetically homozygous. In a previous study, we found that the macropterous and brachypterous strains in three species of rice planthoppers, N. lugens, Sogatella furcifera, and Laodelphax striatellus, did not response in their wing morphs to photoperiod changes [14]. The pure-bred strains of wing morphs showed different responses to environmental changes, as compared to the natural populations [3,14,18]. These results further support the two-stage genetic mode in the wing dimorphism of BPHs ( Figure 9).

Conclusions
We established nearly pure-bred strains of the long-winged and short-winged morph in BPHs and then analyzed the transcriptome of them. Wing morph selection induced significant changes in gene expression levels between the two strains across three developmental stages (egg, nymph, and adult), and gene expression profiles were different. The significantly enriched pathways for DEGs are either constant across three development stages or stage-dependent. We hypothesized that these constant DEGs might be involved in the solely genetic determination of wing morphs, whereas stage-dependent DEGs might be involved in the wing morph determination by gene and environment interactions. Genes expressed stage-independently and stage-dependently may modulate the wing dimorphism of BPHs.