Genome sequence of the small brown planthopper, Laodelphax striatellus

Sample and sequencing

The inbreeding line used for genome sequencing is an inbred laboratory strain that was derived from a field population collected in Hai’an, Jiangsu province, China. A single gravid female was selected, and her progenies were sib-mated for 22 generations to obtain the inbreeding line. Planthoppers were reared on 2–3-cm rice seedlings at 25°C and a photoperiod of 16:8 hours of light/dark. DNA was extracted by using Puregene Core Kit A (Qiagen Sciences, Maryland, USA) from the F22 specimens following the manufacturer's instruction. We built 5 libraries with insert size between 180 bp and 800 bp for paired-end sequencing and 9 libraries with insert size between 1.4 Kb and 24 Kb for mate-pair sequencing according to the standard protocols of the Illumina HiSeq 2500 sequencer (Table 1). We also constructed 33 Pacbio RSII libraries according to the standard Pacbio protocols (Table 1). In total, we generated 190 Gb of Illumina data (126 Gb of paired-end reads and 64 Gb of mated-pair reads) and 32.4 Gb of Pacbio data, representing 316× and 54× coverage of the genome, respectively.

For transcriptome sequencing, total RNA was isolated from 4 tissues (antenna, brain, fatty body, and gonad) and whole bodies of 3 developmental stages (egg, nymph, adult) of L. striatellus using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. Nanodrop (Thermo Scientific, Wilmington, DE, USA) was used to determine RNA quantity, and gel electrophoresis was used to examine RNA quality. cDNA libraries were constructed according to the manufacturer's instructions and sequenced on an Illumina HiSeq 2500 sequencer.

Estimation of genome size and determination of chromosome number

We estimated the genome size of L. striatellus using 2 independent approaches: flow cytometry [10] and k-mer analyses [11]. The flow cytometry analysis was carried out according to a published procedure [10]. Briefly, a female adult was ground in PBS-T buffer. The mixture was filtered by a 40-μm cell filter, incubated with 2 μg/mL RNase A at 37°C for 15 minutes, and then stained with 5 μg/mL propidium iodide at 25°C for 30 minutes. The fluorescence signal was detected by a FACSCallbur Analyzer (Becton-Dickinson, San Jose, CA, USA). Heads of Drosophila melanogaster and cytoblasts of Gallus gallus were treated with the same procedure as genome size references. The genome sizes of D. melanogaster and G. gallus are known to be 0.18 pg and 1.25 pg, respectively [12]. As shown in Fig. S1, the genome size of L. striatellus was estimated to be 0.60 pg (587 Mb) by the flow cytometry method. In k-mer analysis, 31.94 Gb of clean reads were utilized to generate a k-mer (k = 17) depth distribution curve (Fig. S1D), based on which the genome size was estimated to be 550 Mb. Accordingly, the haploid genome size of L. striatellus was estimated to be 550–587 Mb.

The chromosome number was determined by cytological analysis of testes cells. The testes of newly emerged males were dissected in insect Ringer solution, fixed in Carnoy's fixative for 15 minutes. The testes were washed with 0.01 mol/L PBS solution, stained at 0.5 μg/mL Hoechst 33 258, and sealed with Antifade Mounting Medium (Beyotime, Jiangsu, China). Cells in meiosis phase were selected for chromosome counting under a confocal microscope Zeiss LSM710 (Carl Zeiss, Oberkochen, Germany). In most cases, 15 haploid chromosomes were observed (30 for diploid chromosomes) (Fig. S2), although sometimes only 14 were visible. Thus the number of chromosomes in L. striatellus was determined to be 2n = 30.

Genome assembly and assessment

We assembled the genome with both Illumina sequencing and Pacbio sequencing data. Illumina data were used to build contigs and scaffolds as follows. First, all reads with ≥10% unidentified nucleotides, or with >10 nt aligned to the adapter sequences, or being putative PCR duplicates were removed to obtain clean reads. Mate-pair reads from libraries with insert sizes >2 kb were classified as paired-end, unpaired, negative, and mate-pair reads, and only the negative and mate-pair reads were retained for the assembly. Second, we employed SOAPdenovo v. 3.0 (SOAPdenovo, RRID:SCR_010752) [13,14] with the parameters “pregraph -K 33 -p 30 -d 30; contig –k 33 –M 3” to build de Bruijn graph and assemble sequencing reads into contigs. Third, all mate-pair reads were mapped to the contigs, and mate-pair information was added in a stepwise manner to connect contigs into scaffolds. GapCloser v. 1.12 (GapCloser, RRID:SCR_015026) [13] was used to fill the gaps between scaffolds with a local assembly strategy. Afterwards, PBJelly v. 15.8.24 (PBJelly, RRID:SCR_012091) [15] was used to fill the gaps between scaffolds using the 32.4 Gb (∼54×) of Pacbio data. Briefly, all the gaps (length > 25 bp) on the assembly were identified first, and the Pacbio reads were mapped to the assembly using PBJelly. The BLASR alignments were parsed to identify gap-supporting reads by comparing aligned and unaligned base positions within each read [16]. Overlap-layout-consensus engine ALLORA within PBSuite (v. 15.8.24, Pacific Biosciences Menlo Park) [17] was used to assemble the reads for each gap to generate consensus gap-filling sequences. As the final step, the consensus gap-filling sequences were spliced into the corresponding gap position in the draft assembly, replacing all Ns if the gap was closed and leaving the appropriate number of Ns if the gap was only reduced.

With the above assembly procedure, we obtained a final assembly of 541 Mb, having 38 193 scaffolds with a contig N50 length of 118 Kb and a scaffold N50 length of 1.1 Mb. The length of the assembly accounts for 91.7% and 98.4% of the estimated genome size by flow cytometry and k-mer analysis, respectively. The longest contig and scaffold were 2.0 Mb and 10.4 Mb, respectively (Table 2). The Pacbio sequencing data greatly improved the length of contigs compared with the published genomes of N. lugens (contig N50, 24.2 Kb [8]) and S. furcifera (contig N50, 70.7 Kb [9]), which were assembled with Illumina data only (Table 2). We aligned clean reads onto the genome assembly using BWA (BWA, RRID:SCR_010910) [18] and calculated the fraction of bases at a given sequencing depth. The results showed a very small fraction of low-coverage bases, suggesting high coverage and accuracy of the genome assembly (Fig. S3).

Validation and quality control

The completeness and accuracy of the genome assembly were assessed by 4 independent approaches. First, the overall base composition and the percentage of Ns were calculated. As shown in Table S1, the assembled genome had a low percentage (1.99%) of Ns and an expected base composition, which is similar to that of the other 2 planthoppers. The overall GC content of L. striatellus was 34.54%, similar to that of N. lugens [8] and slightly higher than that of S. furcifera [9]. Second, we remapped Illumina paired-end reads to the assembly using BWA [18], and we found that 93.2% of reads could be mapped back, covering 96.83% of the assembled genome, including 95.08% of the genome with ≥×20 coverage (Table S2). Third, we performed de novo transcriptome assembly using Trinity v. 2.0.2 (Trinity, RRID:SCR_013048) for RNA-seq data from multiple developmental stages and tissues (Table 1). We also included 2 published RNA sequencing datasets from salivary glands and alimentary canal [19] in the transcriptome assembly. We mapped the assembled transcripts to the genome assembly using TopHat (TopHat, RRID:SCR_013035) with default parameters and found that 90.31% of the transcripts with >90% transcript coverage were aligned to 1 scaffold (Table S3), indicating that most expressed genes were correctly assembled in the genome. When the RNA reads from the 9 transcriptome datasets were directly mapped to the genome, 78% to 94% could be correctly mapped to the genome with appropriate splicing, indicating that the genome assembly had a good representative of gene regions (Table S4). Finally, the benchmarking universal single-copy orthologs v. 1 (BUSCO, RRID:SCR_015008) dataset representing 2675 genes for arthropoda was used for genome assessment [20]. Our assembled genome captured 92% (2470/2675) of the BUSCO genes, suggesting that a gene repertoire was nearly complete (Table S5). Taken together, these results suggest that our assembled genome was highly accurate and nearly covered the whole genome.

Annotation of repetitive elements

Two independent methods, namely homology-based and de novo prediction, were applied for repetitive element annotation. For the homology-based method, the assembled genome was compared with Repbase, issued on 13 January 2014 [21], using RepeatMasker v. 4.0.5 (RepeatMasker, RRID:SCR_012954) and RepeatProteinMasker (v. 1.36) with default settings [22]. For the de novo prediction, we built a de novo repeat library with LTR_FINDER v. 1.0.5 (LTR_Finder, RRID:SCR_015247) [23], Piler (v. 1.06) [24], RepeatScout v. 1.0.5 (RepeatScout, RRID:SCR_014653) [25], and RepeatModeler v. 1.0.8 (RepeatModeler, RRID:SCR_015027). Tandem Repeat Finder (v. 4.07b) [26] was used to search tandem repeats. Furthermore, RepeatProteinMask [22] was used to identify putative transposable element (TE)–related proteins. After merging all the repetitive elements identified by abovementioned tools, we identified a total of 139.1 Mb of repetitive sequences, accounting for 25.7% of the genome (Table S6). The percentage of repetitive elements in the L. striatellus genome was much lower than those of N. lugens (48.6%) [8] and S. furcifera (44.3%) [9]. Of all the repetitive sequences, 10.59% were the class I transposable elements (retrotransposon), including 5.01% long interspersed nuclear elements, 1.32% long terminal repeats, and 4.26% short interspersed nuclear elements. Class II elements (DNA transposons) represented only 4.92% of the genome (Table 3). L. striatellus had the lowest TE fraction and the smallest genome size compared with N. lugens and S. furcifera (Table 3).

Annotation of protein-coding genes

The protein-coding genes were annotated with evidence from the homology-base method, ab initio prediction, and RNA-seq data. For the homology-based method, the annotated gene sets from 8 species, N. lugens, Acyrthosiphon pisum, Pediculus humanus, Nasonia vitripennis, D. melanogaster, Bombyx mori, Rhodnius prolixus, and Daphnia pulex (Table S7), were aligned to the L. striatellus genome using TBLASTN (TBLASTN, RRID:SCR_011822) [27] with an E-value cutoff of 1E⁻⁵. GeneWise v. 2.2.0 (GeneWise, RRID:SCR_015054) [28] was used to define gene models. For ab initio prediction, we utilized Augustus v. 3.1 (Augustus: Gene Prediction, RRID:SCR_008417) [29], GlimmerHMM v. 3.0.4 (GlimmerHMM, RRID:SCR_002654 [30], SNAP (v. 2013–11-29) [31], GeneID (v. 1.4) [32, 33], and GENSCAN v. 1.0 (GENSCAN, RRID:SCR_012902) [34] to predict potential protein-coding genes from the repeat-masked genome. Furthermore, we identified gene structures with the assistance of 9 transcriptomes assembled by Tophat-Cufflinks (v. 2.2.1) [35] and Trinity-PASA (v. 2.0.2) [36], respectively. Then we integrated all predicted gene structures above with EvidenceModeler (v. 1.1.1) [37] to obtain a nonredundant set of 17 736 protein-coding genes with an average gene length of around 16.17 Kb (Table S8–S9, Fig. S4). We constructed the orthologous gene families using annotated genes from 22 closely related species (Table S7) and found that L. striatellus had 4210 species-specific genes, fewer than those of N. lugens (10 163) and S. furcifera (7743) (Fig. 2). This may be attributed to the smaller genome size and lower gene number in L. striatellus.

We used 3 methods to evaluate the gene models that we obtained. First, we examined the 2 Kb upstream and downstream regions of annotated genes and found that the majority (16 525, 93.17%) of genes did not contain any ambiguous bases (Ns) in the 2 Kb up- and downstream regions, indicating that these gene models are not located near an assembly gap and thus the gene models are unlikely to be a fragment. Second, we compared our annotated genes with the corresponding orthlogous genes in D. melanogaster. We performed BLASTX (BLASTX, RRID:SCR_001653) [27] searches against the D. melanogaster gene set using the de novo assembled transcripts in L. striatellus. A total of 8484 assembled transcripts that had identity >60% with a D. melanogaster gene and covered >90% of the coding region were regarded as full-length transcripts. Among them, 3728 transcripts (excluding redundant protein isoforms) containing a complete ORF were searched against the annotated genes, and 3093 (82.97%) of them had a near perfect match to an annotated gene, indicating that most annotated genes were complete. Third, we compared our annotated genes to the 2 sets of BUSCO (v. 2) genes (1066 arthropoda genes and 303 eukaryote genes) [20] and found that our predicted genes were considered complete BUSCO genes in 97.6% and 98.0% of the eukaryote genes and arthropoda genes, respectively (Fig. S5), suggesting that a nearly complete repertoire of protein-coding gene set was determined.

To estimate the level of heterozygosity in the gene model, we aligned 23× reads to the genome assembly with BWA [18]. After removing duplicates, heterozygous SNPs were identified using BCFtools [38]. The heterozygous SNPs in the coding regions of each gene were used to compute read coverage and heterozygosity. Only a single heterozygosity peak of around 0.3 was detected (Fig. S6A). We ranked the heterozygosity rate of all the gene set and took the top 20% as high heterozygosity (the remainder was designated low heterozygosity). Coverage histograms of high and low heterozygosity showed similar ranges of coverage distribution (Fig. S6B). Therefore, the heterozygosity did not influence the gene annotation.

In order to obtain putative functional assignments to the annotated genes, we compared the annotated protein sequences of L. striatellus to proteins in the Kyoto Encyclopedia of Genes and Genomes (KEGG, RRID:SCR_012773) [39], NR [40], and Swiss-Prot [41] databases using BLASTP (BLASTP, RRID:SCR_001010) [27] with an E-value cutoff of 1E⁻⁵. Domains and motifs were scanned in Interpro [42] database by InterProScan (InterProScan, RRID:SCR_005829) [43]. There were 78.7%, 66.3%, 63.6%, and 69.5% of annotated proteins showing significant sequence similarity with the proteins in NR, Swiss-prot, KEGG, and InterPro (InterPro, RRID:SCR_006695), respectively. Among the 12 322 genes with an InterPro hit, 11 159 (90.6%) had Pfam (Pfam, RRID:SCR_004726) annotations and 8935 (72.5%) had gene ontology (GO, RRID:SCR_002811) associations. After removing redundancy, 14 182 of 17 736 genes (80.0%) were assigned to known databases (Fig. 3). Among the 3554 unannotated genes, 1391 (7.8%) were L. striatellus–specific genes.

Gene orthology prediction

Twenty-one sequenced insects (Zootermopsis nevadensis, Tribolium castaneum, Anoplophora glabripennis, Anopheles gambiae, D. melanogaster, A. pisum, Diuraphis noxia, Cimex lectularius, L. striatellus, R. prolixus, N. lugens, S. furcifera, Diaphorina citri, Oncopeltus fasciatus, Apis mellifera, N. vitripennis, B. mori, B. tabaci, Danaus plexippus, Locusta migratoria, and P. humanus) and 1 noninsect arthropoda sequenced species (D. pulex) were used to infer gene orthology and reconstruct the phylogenetic tree. The annotated coding sequences were downloaded from the websites listed in Table S7. The homologous gene families were identified using TreeFam [44, 45] and ascribed in different categories (Fig. 2). The gene families were identified following these steps: (i) BLASTP [27] was used to compare all protein sequences for the 22 species with an E-value cutoff of 1E⁻⁷; (ii) the blast alignments were concatenated by Solar (v. 0.9.6) [45], followed by homology identification among protein sequences; and (iii) gene families were identified using hcluster_sg (v. 0.5.0) [45]. RAxML (v. 8.0.19) [46] was used to reconstruct the phylogenetic tree based on the concatenated single-copy protein sequences under the PROTGAMMAAUTO model with 100 bootstrap replicates. R8s (v. 1.7.1) [47] and MCMCtree (PAML package, v. 4.7; PAML, RRID:SCR_014932) [48] were used to estimate the divergence times among species. The parameters used in MCMCtree were “–rootage 510 -clock 3 -alpha 0.977999 -model 7.” To examine gene family expansion and contraction in the 3 planthoppers, we chose 1 additional hemipteran species, R. prolixus, as outgroup to infer expanded/contracted gene families using CAFE (v. 3.1) [49]. A conditional P-value was calculated for each gene family, and the gene families with P-values <0.05 were considered as significantly expanded or contracted. The phylogenetic analysis revealed that L. striatellus clustered together with the other 2 planthoppers and had a closer relationship to S. furcifera than N. lugens (Fig. 4). The divergence times of nonplanthopper insect species were generally consistent with those estimated in the previous study [8]. The results of molecular dating analysis indicated that the ancestor of L. striatellus and S. furcifera split with N. lugens about 87.5 million years ago and that L. striatellus diverged from S. furcifera approximately 31 million years ago (Fig. S7).

Compared with N. lugens and S. furcifera, L. striatellus had fewer expanded gene families and more contracted gene families (Fig. S8). This might partially explain why L. striatellus has the lowest gene number among the 3 planthopper species. Since the divergence of L. striatellus and S. furcifera, L. striatellus and S. furcifera have had 95 and 547 expanded gene families, respectively (Fig. S8). The significantly expanded gene families in L. striatellus included some specific members of multigene families, such as odorant receptor, cytochrome P450, and serine protease (especially trypsin) (Table S10). The specific members of chemosensory protein, odorant binding protein, carboxylesterase, and ATP-binding cassette transporter families were also increased in L. striatellus although their P-values were higher than 0.05 (Table S10). Expansion of these gene families may have contributed to the widest host plant range of L. striatellus among the 3 planthoppers. The specific members of gene families associated with energy metabolism were significantly expanded in S. furcifera, such as acyl-CoA synthetase, fatty acyl-CoA reductase, acyl-CoA-binding protein, and acyl-coenzyme A thioesterase. The specific members of glyceraldehyde-3-phosphate dehydrogenase, D-beta-hydroxybutyrate dehydrogenase, ADP/ATP translocase, acyl-CoA transporter, and ATP synthase families also increased, although with P-values higher than 0.05 (Table S10). N. lugens had 433 expanded gene families (Fig. S8). A bunch of specific members from energy metabolism–related gene families, including Delta(3,5)-Delta(2,4)-dienoyl-CoA isomerase, ATP-citrate synthase, malonyl-CoA decarboxylase, NADH dehydrogenase (ubiquinone) 1α subcomplex subunit 7 and subunit 8, acyl-CoA synthetase, ATP synthase, and enoyl-CoA delta isomerase increased in N. lugens although their P-values were higher than 0.05 (Table S10). Expansion in the energy metabolism–related gene families is in accordance with the migratory habits of S. furcifera and N. lugens.

Olfaction and detoxification system

It is essential for herbivorous insects to recognize and locate their host plants utilizing their sense of gustation and olfaction. Chemicals from the environment are received and recognized by chemoreceptor genes, including odorant receptors (ORs), gustatory receptors (GRs), and ionotropic receptors (IRs) in gustatory and olfactory organs. Detoxification gene families also play an essential role in defense against natural xenobiotics from host plants or synthetic xenobiotics including insecticides. To identify chemoreception and detoxification-related genes in L. striatellus, we retrieved corresponding gene sequences of other insect species from previous studies and used them as queries. These genes were searched against the L. striatellus gene set using BLASTP [27] with an E-value cutoff of 1E⁻⁵. In addition, we scanned the gene sets of 3 planthoppers for domain information using InterProScan and extracted genes with domains corresponding to each family. Finally, we integrated results from both BLASTP and InterProScan to obtain the final set of protein families.

There were 106 ORs, 38 IRs, and 12 GRs identified in L. striatellus (Table S11). The numbers of ORs and GRs in L. striatellus were more than twice as many as those in N. lugens and S. furcifera, representing a significant expansion in these 2 families. This is consistent with the fact that L. striatellus is the most polyphagous among the 3 planthoppers because polyphagous insects tend to have more OR genes than monophagous [8]. Moreover, we identified 2 protein families important for odor recognition and pheromone perception, namely odorant binding proteins (OBPs) and chemosensory proteins (CSPs). There were 16 OBPs and 31 CSPs in L. striatellus, the most among the 3 planthoppers (Table S11). The relatively higher number of odor-related genes in L. striatellus might be closely related to its polyphagous habit.

We manually annotated families of detoxification-related genes, including 26 UDP-glycosyltransferases, 29 glutathione-S-transferases, 54 carboxyl/cholinesterase, 73 ATP-binding cassette transporters, and 76 cytochrome P450s in L. striatellus (Table S12). The total number of detoxification-related genes in L. striatellus was smaller than that in N. lugens, but larger than that in S. furcifera.

Immune-related genes

We identified immune gene repertoires of the 3 planthoppers using a homology-based method. Immune genes from D. melanogaster, A. gambiae, Aedes aegypti, and Culex quinquefasciatus were downloaded from ImmunoDB [50]. Gene sets from the 3 planthoppers were used as queries and searched against the immune genes of the 4 insects, respectively, using BLASTX with an E-value cutoff of 1E⁻⁵. The best hits were selected for further domain architecture analysis using InterProScan and then were confirmed manually. The number of immune-related genes in L. striatellus was 330, which was more than that in N. lugens (289) and S. furcifera (280) (Table S13). The redundant copies of immune genes in L. striatellus mainly included autophagy genes, 1,3-beta-D glucan binding protein genes, clip-domain serine protease genes, and genes of small RNA regulatory pathway members. However, the numbers of C-type lectin genes and Toll-like receptor genes were lower in L. striatellus compared with the other 2 planthoppers.

Transcriptomic responses of 3 planthoppers to their borne plant viruses

L. striatellus, S. furcifera, and N. lugens transmit different rice viruses. To explore the molecular response to respective plant viruses, we analyzed and compared the transcriptomic responses of L. striatellus to RSV, S. furcifera to SRBSDV, and N. lugens to RRSV. The 3 viruses are transmitted in a persistent-propagative way. For L. striatellus, RSV was incubated in the fourth-instar nymphs for 5 days, as described previously [51]. Three replicates of infected or noninfected insects were used to construct paired-end RNA-seq libraries for sequencing on an Illumina Hiseq 2500 sequencer. The transcriptomic data of S. furcifera infected with SRBSDV were retrieved from a previous study [52]. The third-instar nymphs of N. lugens were infected by RRSV for 7 days before being collected for RNA extraction using the SV Total RNA Isolation System (Promega, Madison, WI, USA). The gene expression libraries for RRSV-infected and noninfected samples were constructed and sequenced on an Illumina HiSeq 2000 sequencer. RNA-seq reads were mapped to the corresponding genome using TopHat2 (v. 2.1.1) [53]. For L. striatellus and S. furcifera, HTSeq [54] was used to count the number of reads mapped to each gene model, and the edgeR package was used to identify differentially expressed genes (DEGs) with a fold change cutoff of 2 and FDR cutoff of 0.01. For N. lugens, generalized fold change for ranking differentially expressed genes from RNA-seq data was used to detect DEGs without biological replicates. The gene annotation files were downloaded from the corresponding websites (Table S14). We referred to genes with higher expressions in the viruliferous group as upregulated genes and lower as downregulated. The results showed that 460 (185 up and 275 down), 162 (48 up and 114 down), and 1070 (515 up and 555 down) genes were differentially expressed in L. striatellus, N. lugens, and S. furcifera, respectively, when bearing their respective plant viruses.

The DEGs in the 3 planthoppers were compared in GO terms, and the common GO terms were retrieved (Table S15). The upregulated genes in the 3 planthoppers were involved in the biological processes of regulation of transcription (GO:0 006355) and protein phosphorylation (GO:0 006468). The downregulated genes in the 3 planthoppers took part in the biological processes of carbohydrate metabolic process (GO:0 005975), chitin catabolic process (GO:0 006032), and proteolysis (GO:0 006508).

Two zinc finger proteins of L. striatellus, 1 zinc finger protein of N. lugens, and 6 zinc finger proteins of S. furcifera were commonly upregulated while genes of chitinases, cytochrome P450 CYP4s, and trypsins were commonly downregulated in the 3 planthoppers (Table S16) in response to their respective plant viruses. We also identified homologous genes that were commonly regulated in the 3 planthoppers by aligning N. lugens DEGs with those of L. striatellus and S. furcifera using BLASTP with a cutoff of 1E⁻³, a sequence identity higher than 60%, and a coverage higher than 50%. Three groups of homologous genes, including 1 group of commonly upregulated genes and 2 groups of commonly downregulated genes, were retrieved from the 3 planthoppers (Table S17). The protein lengths of these homologous genes ranged from 120 to 472 amino acids. We used these proteins as queries to search the NR database and found no homologous genes in other species with a cutoff of 1E⁻⁷, indicating that these genes are likely planthopper-specific genes.

Differences in immune response to virus infection in the 3 planthoppers were also observed. The RNAi pathway genes, RISC-loading complex TARBP2 and argonaute-3, were upregulated in S. furcifera and N. lugens, respectively, but genes in the RNAi pathway did not respond to virus infection in L. striatellus. The antimicrobial peptide defensin was upregulated in L. striatellus and N. lugens but was downregulated in S. furcifera. The expression of the Down syndrome cell adhesion molecule gene increased in L. striatellus [55] and decreased in S. furcifera, but did not show significant change in N. lugens in response to their respective plant viruses.

In summary, we reported a high-quality of genome of L. striatellus, a notorious rice pest insect. L. striatellus has the smallest genome and the lowest number of protein-coding genes compared with the other 2 rice planthoppers, S. furcifera and N. lugens. Comparative genomic analyses identified expansions and contractions in olfactory genes, detoxification genes, immune genes, and energy metabolism genes among the 3 rice planthoppers, which may have contributed to their differences in important traits such as host range, migratory habit, and plant virus transmission. Despite having the smallest genome, L. striatellus has the widest host plant range among the 3 planthoppers. This situation is different from that of the genome evolution in Aphididae, where the soybean aphid, Aphis glycines, which is an extreme specialist, has the smallest genome compared with another 3 aphid species with published genome sequences [56]. With the addition of the L. striatellus genome, the genome data of the 3 rice planthoppers will aid studies in various areas of planthoppers and promote control strategies in the future.