\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Chromosome-level genome assembly and population genomic analysis provide
novel insights into the immunity and evolution of \emph{Sogatella
furcifera}}



\author[1]{Ting Cui}%
\author[1]{Qing Bai}%
\author[1]{Weiqiang Yu}%
\author[1]{Dongyang Guo}%
\author[1]{Yawen Ban}%
\author[1]{Kun Chen}%
\author[1]{Ali  Raza}%
\author[2]{Guohui Zhou}%
\author[1]{Qingfa Wu}%
\affil[1]{University of Science and Technology of China}%
\affil[2]{South China Agricultural University}%


\vspace{-1em}



  
  \date{May 21, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
\emph{Sogatella furcifera} is an agricultural pest of great concern in
China and Southeast Asian countries. However, the lack of accurate and
complete reference genome resources has hindered the understanding of
immunity and evolution of \emph{S. furcifera}. Here, we utilized
Nanopore sequencing to generate a chromosome-level assembly and
annotation of the \emph{S. furcifera} genome (0.64 Gb), with a GC
content of 34.25\%. This genome comprised 15 chromosomes covering
95.04\% of the estimated genome size, together with an additional 624
small scaffolds making up the remaining 4.96\% of the genome of \emph{S.
furcifera}. A total of 24,669 protein-coding genes as well as 1211 long
noncoding RNA and 7595 circular RNA transcripts were well annotated and
predicted. Comparative genomic analysis revealed the rapidly evolved
genes associated with multiple immune-related pathways in \emph{S.
furcifera}, which may be responsible for its rapid evolutionary
adaptation. In addition, 14 immune genes in the classic immune pathways
were selected for functional validation through RNA interference
experiments, demonstrating the antiviral effects of Dorsal and Dif genes
in \emph{S. furcifera}. Genome resequencing of 44 individuals from 12
geographic populations revealed an absence of population structures and
frequent gene flow among all populations. Sweep analysis indicated that
2926 genes were under natural selection and significantly enriched in
several biological processes of morphogenesis and immunity. The first
systematic identification of immune genes and noncoding RNAs from
chromosome-level genome assembly plus the comparative and population
genomic analysis will provide more insights into the understanding of
the immunity and evolutionary adaptation of \emph{S. furcifera}.%
\end{abstract}%



\sloppy


\textbf{INTRODUCTION}

The white-backed planthopper (WBPH),\emph{Sogatella furcifera} , can
either directly threaten the rice plant through phloem sap-sucking or
indirectly transmit devastating rice viruses, such as southern rice
black-streaked dwarf virus (SRBSDV), triggering a double blow for rice
production. The WBPH transmits SRBSDV in a persistent circulative and
propagative manner, resulting in one of the largest threats to rice
production in southern China and Southeast Asian countries in the past
decade (Hoang et al., 2011; Matsukura et al., 2013; Zhou, 2010; Zhou,
Xu, Xu, \& Zhang, 2013). Yield losses can reach to 80\% or even more in
local rice fields seriously affected by WBPH and SRBSDV (Liu, 2010).
Therefore, it is necessary to deeply understand the mechanisms of the
WBPH associated with transmitting the virus SRBSDV.

\emph{Sogatella furcifera} is a long-distance migratory pest that cannot
tolerate cold weather and thus migrates long distances northward or
northeastward, following a northeasterly flowing air current each year
from tropical and subtropical regions, including southern China and
Vietnam (Turner, Song, \& Uhm, 1999). \emph{S. furcifera} populations
migrating from the Indo-China Peninsula are considered the primary
source of \emph{S. furcifera} in northern China (Shen, 2003),
whereas\emph{S. furcifera} populations migrating from northern Vietnam
or southern China are considered the primary sources of \emph{S.
furcifera}in northeastern China, Korea and Japan (Matsumoto et al.,
2013). Although the migration, genetic diversity and structure of
\emph{S. furcifera} populations have been investigated across some Asian
countries, the insufficient number and robustness of genetic markers,
including microsatellites, mitochondrial genes, and SNP markers, still
have drawbacks for understanding population features (Liu, Gui, \& Li,
2010; Mun et al., 1999; Nam, Kim, \& Lee, 2019; Sun, Jiang, Wang, \&
Hong, 2014; Xie et al., 2014; Yang et al., 2020). In recent years,
pioneering genomic technologies have provided more resources to screen
highly representative markers for population studies, such as the
genome-wide identification of SNP markers. However, far less information
has been gained from a genome-wide study of the WBPH populations due to
the lack of a high-quality WBPH genome sequence.

Since 2017, the first reference genome of WBPH has been published and
systematically elucidated the cellular pathways and molecular regulatory
networks involved in functional regulation in the WBPH during each
developmental stage (Wang et al., 2017; Wang et al., 2016). Given that
the early data and analyses were based on next-generation sequencing
platforms and the incompleteness of the genomic and transcriptomic
information due to the short fragmentation of read lengths, more
accurate annotation of functional genes and noncoding RNAs is needed to
understand the immune mechanisms of WBPH. Later, another research group
concatenated the first version of the \emph{S. furcifera} genome to
generate a chromosome-level genome by Hi-C scaffolding technology, and
they further updated this genomic data using the Nanopore sequencing
platform (ASM1714138v1) (Ma et al., 2021). All these previous genomic
data of \emph{S. furcifera} have provided important insights to study
the pest disease of rice associated with SRBSDV and WBPH\emph{.}However,
these genomic data were still mis-spliced to some extent, and the
protein-coding genes and noncoding RNAs were not well annotated. This
has hindered the deep and effective study of the immune interaction
between WBPH and SRBSDV, which has been well studied in holometabolous
insect species, especially in dipterans (e.g., fruit flies and
mosquitoes) (Ferreira \selectlanguage{ngerman}Á et al., 2014; Gao et al., 2021; Nazzi et al.,
2012; Xi, Ramirez, \& Dimopoulos, 2008; Zambon, Nandakumar, Vakharia, \&
Wu, 2005; Zhang et al., 2020). Therefore, there is a great need to
provide more accurate and complete information for a deeper
understanding of the immune mechanisms of the WBPH.

In this study, two mainstream third-generation sequencing platforms were
applied to obtain more accurate and complete genomic and transcriptomic
data of WBPH. Here, we reported the high-quality genome and
transcriptome sequences of \emph{S. furcifera} , adding this important
species to the other sequenced insects of the order \emph{Hemiptera} .
In addition, the resequencing of 44 individuals collected from 12
geographic locations in Vietnam and China revealed the genome-wide
variation of \emph{S. furcifera} populations. Taken together, these
information provide more genomic insights into understanding of the
immunity and population evolution of \emph{S. furcifera} .

\textbf{METHODS}

\textbf{Samples, DNA and RNA extraction}

The white-backed planthopper sequenced in the present study originated
from the University of Science and Technology of China. A laboratory
colony was maintained at 26°C with 70\% humidity under a 16:8 h
light/dark photoperiod for two years, spanning at least 20 generations.
An inbreeding line was obtained by single pair sib mating for 18
generations. To assist gene annotation, we performed whole transcriptome
sequencing. Total RNA was prepared from \emph{S. furcifera} specimens at
all developmental stages of the white-backed planthopper (embryos,
1\textsuperscript{st}--5\textsuperscript{th}instar nymphs, 5-day-old
adults and 10-day-old adults) using TRIzol reagent (Invitrogen). RNA
quantitation was performed by UV absorbance, and its quality was further
confirmed by gel electrophoresis. The 44\emph{S. furcifera} individuals
were collected across China and Vietnam from 12 geographic populations
in 2016.

\textbf{Library construction and high-throughput sequencing}

The genomic library for Oxford Nanopore Technologies (ONT) sequencing
was constructed using the SQK-LSK109 ligation sequencing kit according
to the manufacturer's instructions, and the libraries were subsequently
sequenced using an Oxford Nanopore PromethION sequencer by Biomarker
Technologies Corporation (Beijing, China). The RNA library was prepared
using the SMARTer\selectlanguage{english} PCR cDNA Synthesis Kit and sequenced on the PacBio
sequencing platform by Biomarker Technologies Corporation (Beijing,
China).

\textbf{Contig-level genome assembly}

The genome was primarily assembled based on the Nanopore sequencing data
(\textbf{Table S1} ), followed by three round corrections of the genome
assembly using Racon software
(https://hpc.ilri.cgiar.org/racon-software) and then one round
correction of the genome assembly using Medaka software
(https://hpc.ilri.cgiar.org/medaka-software). The process was continued
with three rounds of correction based on Illumina reads using Racon
software, with each round being evaluated by benchmarking universal
single-copy orthologs (BUSCO) to determine if it can be used as a contig
level sketch before scaffold linkage. Another software (NextPolish) was
further used for genome correction (Hu, Fan, Sun, \& Liu, 2019). For
this data set, the NextPolish correction process eliminated fewer
sequences, with a higher corrected genome size and better BUSCO score
than those of Racon. The correction step greatly facilitated the
completeness of the contig-level genome assemblies through the
NextDenovo assembly combined with NextPolish correction based on
Nanopore reads as well as NextPolish based on Illumina reads.

\textbf{Chromosomal-level genome assembly using Hi-C data}

As \emph{S. furcifera} is sex-determined by the XO chromosome system and
females yield nearly twofold larger sequencing coverage than males, the
female adult WBPH was collected in this study for Hi-C library
construction following a standard procedure (Xie et al., 2015). In
practice, Hi-C was performed in two rounds for library construction and
sequencing, yielding a total of 264,745,379 reads. The Hi-C reads were
mapped to the contig level genome with Juicer software (Durand et al.,
2016), and duplicate information records were obtained after the
comparison, filtering and deduplication processes. The contig level
genome was then concatenated to the scaffold level by 3D-DNA software
(Dudchenko et al., 2017) based on this matching information. As the
genome file contains inversions, translocations and unclear chromosome
boundaries, a manual review and refinement step for quality control and
interactive correction were conducted by loading Hi-C interactions in
the Juicebox Assembly Tools (JBAT) (Dudchenko et al., 2018). A
correlation heatmap was presented to illustrate the interaction of each
chromosome, and then a manual correction of the errors was performed to
obtain a corrected position coordinate file. Finally, the genome of the
intermediate steps was adjusted based on this coordinate file to obtain
the final polished chromosomal-level reference genome. In total, 15
chromosome groups for \emph{S. furcifera} assembly were obtained.

\textbf{Genome annotation}

Protein-coding gene annotation

The structural and functional genome annotation followed pipeline B of
the BRAKE software (Hoff, Lange, Lomsadze, Borodovsky, \& Stanke, 2015),
combined with the genomic file and transcriptomic data. The raw reads
from RNA sequencing of all eight different developmental stages were
mapped onto the \emph{S. furcifera} genome assembly, and alternative
spice site information was obtained. These transcripts were used to
retrain GeneMark-ET, which is a self-training algorithm (Lomsadze,
Burns, \& Borodovsky, 2014), and SNAP (Korf, 2004). Thereafter, the same
mapping information was used for gene model-based prediction through
AUGUSTUS software (https://github.com/Gaius-Augustus/BRAKER) (Stanke \&
Morgenstern, 2005). Finally, all predicted gene structures were
integrated to yield a consensus gene set followed by the automated MAKER
pipeline (Cantarel et al., 2008). Annotated genes were used to obtain
gene function, and their protein sequences were aligned against
different protein databases, including eggNOG (Insecta) (Jensen et al.,
2008), PANNZER2 (arthropoda)
(http://ekhidna2.biocenter.helsinki.fi/sanspanz/), Interproscan (Jones
et al., 2014) and KAAS (40 organisms of Insecta) (Moriya, Itoh, Okuda,
Yoshizawa, \& Kanehisa, 2007). Finally, the GO (Gene Ontology), KEGG
(Kyoto Encyclopedia of Genes and Genomes) and Pfam (Protein families
database) annotations and possible pathways of all annotated genes were
obtained.

Noncoding RNA annotation

Noncoding RNAs, including long noncoding RNAs (lncRNAs) and circular
RNAs (circRNAs), were also identified from the assembled transcripts
(\textbf{File S1} ). For the identification of lncRNAs, the pipeline was
showed in \textbf{Figure S1} . For the identification of circRNAs, the
circRNA detection tools `find\_circ' and `CIRI2' were employed in this
study (\textbf{File S1} ). Finally, heatmaps were generated to show the
differential expression of identified lncRNAs and circRNAs at eight
different developmental stages.

DNA repeat annotation

The repeat sequences were also identified through alignment with the
known repeat sequences in the databases (Repbase and Dfam) using the
RepeatMasker program (https://www.repeatmasker.org/) with the AB-Blast
search engine (http://www.advbiocomp.com/blast.html). Subsequently, the
species-specific repeat sequences were identified through the initial
annotation based on the structural features of the repeat sequences
using RepeatModeler2 (Flynn et al., 2020). Finally, all the repeat
sequences obtained from the above steps were combined and marked using
lowercase letters.

\textbf{Comparative genomics}

Identification of orthologous genes and phylogenetic analysis

The orthologous and paralogous gene clusters were identified through
OrthoFinder (Emms \& Kelly, 2019) in the assembled genomes of \emph{S.
furcifera} and other 15 related insect species (\textbf{File S1} ).
Additionally, the common water flea (\emph{Daphnia pulex} ) was chosen
as an outgroup during the inference of the phylogenetic species tree
(\textbf{File S1} ). Prior to the construction of the phylogenetic tree,
the optimal model of protein evolution was selected using ProtTest
v3.4.2 (Darriba, Taboada, Doallo \& Posada, 2011) with the parameters
``-all-distributions -F -AIC -BIC -tc 0.5,'' and the best model was
selected as ``LG+I+G+F'' based on BIC. Subsequently, the
maximum-likelihood (ML) phylogenetic species tree was inferred using
RAxML v8.2.10 (Stamatakis, 2014) with the best fit model
PROT+GAMMA+ILGF, and rapid bootstrap inference was implemented with 500
duplications.

To infer the divergence times of different species in the phylogeny, the
MCMCtree program within the Phylogenetic Analysis by Maximum Likelihood
(PAML) package (Yang, 2007) was employed for estimation of divergence
time with the JC69 model. Additionally, the species divergence time was
calibrated with seven calibration time points: 223.8--344.7 million
years ago (Mya) for \emph{D. melanogaster} and \emph{Danaus plexippus} ,
280.9--361.6 Mya for \emph{D. melanogaster} and \emph{Anoplophora
glabripennis} , 312.9--389.7 Mya for \emph{D. melanogaster} and \emph{A.
mellifera} , 172.6--373.3 Mya for \emph{Pediculus humanus}
and\emph{Bemisia tabaci} , 234.9--366.2 Mya
for\emph{Nilaparvata\_lugens} and \emph{Rhodnius\_prolixus} ,
112.5--391.7 Mya for \emph{Bemisia tabaci} and \emph{Rhodnius prolixus}
, and 172.6--373.3 Mya for \emph{Pediculus humanus} and \emph{Bemisia
tabaci} . The calibration time points were collected from the fossil
records and the TimeTree website (http://www. timet ree.org), and the
parameters of mcmctree were set as follows: burn-in=5000,
sample-number=20000, and sample-frequency=50. i-TOL
(https://itol.embl.de/) was used to visualize the phylogenetic tree.

Gene family expansion and contraction and positively selected gene
analysis

Computational Analysis of gene Family Evolution (CAFE) was used to
detect gene family expansion and contraction (De Bie, Cristianini,
Demuth \& Hahn, 2006), in which gene family change was simulated using a
stochastic birth and death model. The gene count file and an ultrametric
species tree with branch length information were provided as the input
for CAFE analysis. The optimal lambda parameter was automatically
determined by the program. For the CAFE analysis result, the gene
families with a family-wide \emph{P} value \textless{} 0.05 were defined
as rapidly evolving families. Moreover, the site models and likelihood
ratio tests (LRTs) were applied for the positive selection analysis
(\textbf{File S1} ).

Immune gene identification and RNAi effects of representative genes in
classic immune pathways

To better understand the evolution of immune genes, they were further
identified and annotated across the WBPH genome based on a homology
search of WBPH protein-coding genes against immune-related genes of
model species (\emph{Drosophila melanogaster} ) and ImmunoDB
(http://cegg.unige.ch/Insecta/immunodb). Additionally, 14 representative
genes (\emph{Toll} , \emph{Dorsal} and \emph{Dif} , \emph{IMD}
and\emph{Relish} , \emph{Hop} and \emph{STAT} , \emph{Atg13}
and\emph{Vps15} , \emph{Tep1} and \emph{Eater} , \emph{Diap1}
and\emph{Dcp1} ) in RNAi, IMD, Toll, Jak-STAT, autophagy, phagocytosis,
and apoptosis pathways, respectively, were characterized and
functionally validated through RNAi knockdown. Double-stranded RNAs
(dsRNAs) were designed by SnapDragon
(https://www.flyrnai.org/cgi-bin/RNAi\_find\_primers.pl) using coding
sequences of each gene to be knocked down. DNA templates for dsRNA
synthesis were amplified with primers containing the T7 RNA polymerase
promoter at both ends. The dsRNA for GFP was used as a negative control
for the nonspecific effects of dsRNA. 2\textsuperscript{nd-} to
3\textsuperscript{rd}-instar \emph{S. furcifera} nymphs (100 to 200)
were used for each gene treatment, and each treatment was conducted in
triplicate. Each nymph was injected with 20--40 nl of dsRNA at a
concentration of 5 \selectlanguage{greek}μ\selectlanguage{english}g/\selectlanguage{greek}μ\selectlanguage{english}l. After injection, the nymphs were reared in
tissue culture jars with rice seedlings, each of which was sealed with
nylon nets and contained at most 100 nymphs. One day after injection,
samples were collected from a set of 5 nymphs to evaluate the RNAi
effects of each gene using qPCR with designed primers.

\textbf{Population genomics analysis}

Genomic resequencing

The re-sequenced strains of \emph{S. furcifera} , consisting of 44
individuals, were collected across China and Vietnam from 12 geographic
populations of \emph{S. furcifera} , including three populations from
Vietnam, Thai Binh (VNTB), Bac Giang (VNBG), and Nam Dinh (VNND), and
nine other populations from China, including Guangxi (GXYZ),
Guangdongguangzhou (GDGZ), Guangdonghaifeng (GDHF), Yunnanshidian
(YNSD), Yunnanmangshi (YNMS), Yunnanchengjiang (YNCJ), Hainanbaoting
(HNBT), Hainanlingshui (HNLS), and Hainanqionghai (HNQH). Genomic DNA
was extracted using the Blood \& Cell Culture DNA Midi Kit (Qiagen) for
genome sequencing. A total of 111 Gb of data were obtained, and the
sequencing depth of each individual sample reached 4- to 5-fold and
covered approximately 95\% of the genome.

Single-nucleotide polymorphism (SNP) calling

The reads from each re-sequenced accession were aligned to the reference
genome using Burrows-Wheeler Aligner (BWA) (Li \& Durbin, 2009). After
mapping, SNP variants were identified based on GenomeAnalysisTK (McKenna
et al., 2010). Variants of each isolate were merged into a single genome
Variant Call Format (gVCF) file using GATK CombineGVCFs. Then, the
combined file was filtered with parameter settings as ``QD\textless{}2.0
\textbar{}\textbar{} FS\textgreater{}60.0 \textbar{}\textbar{}
MQ\textless{}40.0 \textbar{}\textbar{} MQRankSum\textless{}-12.5
\textbar{}\textbar{} ReadPosRankSum\textless{}-8.0 \textbar{}\textbar{}
SOR\textgreater{}3'' according to the recommendation in the GATK Best
Practice Pipeline with modifications
(https://github.com/gatk-workflows/). Joint genotyping was implemented
on the combined file using GATK GenotypeGVCFs. The obtained SNPs were
filtered using the following criteria: genotype missing rate=0 and MAF
(minor allele frequency) \textgreater{} 0.05. Additionally, insertions
and deletions (InDels) were not considered during our downstream
analysis.

Population genetic structure

To ensure that adjacent SNPs were not linked, the squared correlation
(r\textsuperscript{2}) was calculated to estimate the linkage
disequilibrium between each pair of SNPs. An SNP was filtered out when
r\textsuperscript{2} was greater than 0.2. The remaining 465,882
unlinked SNPs were used to construct a phylogenetic tree based on Nei's
genetic distance with 1000 replications using the function `aboot'
within the `poppr' package in the R program 4.2.2. Moreover, the genetic
relationships among individual multilocus genotypes were also analyzed
using the `bitwise.dist' function and visualized by a minimum spanning
network within the `poppr' package in the R program 4.2.2. Additionally,
TreeMix analysis was conducted to reveal possible gene flow events among
the 12 geographic populations (Pickrell \& Pritchard, 2012).

To detect the population structure, discriminant analysis of principal
components (DAPC), a non-model-based method developed and implemented in
the `adegenet' R package, was performed to generate a scatterplot of the
membership probability of each isolate in the corresponding genetic
lineages (Jombart, 2008; Jombart et al., 2010). The `find.clusters'
function in the `adegenet' package was first run on the SNP data set for
different values of K representing the latent population number, and an
\selectlanguage{greek}α \selectlanguage{english}score optimization was used to determine the number of principal
components to retain. A K value is considered good if it exhibits a
lower BIC value than others. To assess how the spatial populations
differed from each other, DAPC was also performed based on the
predefined populations based on different geographic locations. DAPC was
carried out for all SNP loci, and \selectlanguage{greek}α \selectlanguage{english}score optimization was used to
determine the number of principal components to retain.

Genome-wide scan signal of selective sweeps

For sweep analysis, we identified selective sweeps using SweeD (Sweep
Detector) by analyzing the site frequency spectra (Pavlidis et al.,
2013), an open-source tool for the rapid detection of selective sweeps
in whole genomes. The composite likelihood ratio (CLR) test statistics
in SweeD were calculated with default parameters with 10 kb intervals.
High CLR values denote a candidate region for a selective sweep (Nielsen
et al., 2005).

Functional enrichment analysis

The relevant functions of the protein-coding genes overlapping with
selective sweeps were characterized by searching for overrepresented GO
(gene ontology) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes)
pathways (Kanehisa et al., 2016) using the package clusterProfiler.
The\emph{Sogatella furcifera} protein sequences annotated in this study
were used to perform functional enrichment tests on the target genes.
KEGG pathways with a \emph{P} value \textless{} 0.05 were considered
statistically significantly enriched.

\textbf{RESULTS}

\textbf{Genome assembly}

Contig-level genome assembly of Sogatella furcifera

Four software programs, Wtdbg2, NextDenovo, Racon and Nextpolish, were
used to perform the genome assembly and correction (\textbf{Figure 1a,
Table S2} ). The BUSCO of the genome assembled by Wtdbg2 was 77.7\%, and
after multiple rounds of correction, the BUSCO reached 96.3\%, the
repost rate of the Illumina sequencing library reached 98.79\%, and the
final contig level was 716 Mb. The BUSCO of the genome initially
assembled by NextDenovo was 93.0\%, and the genome size was 640.69 Mb.
After multiple rounds of correction, it was found that the BUSCO genome
no longer increased since the first-round correction of the Illumina
sequencing data, generating the final contig level genome with a size of
624.32 Mb. The raw contig-level assembly comprised 443 contigs with a
contig N50 of 4,175,853 base pairs (bp) (\textbf{Table S2} ).

Chromosomal-level genome assembly using Hi-C data

The Hi-C data were further used to improve the genome assembly
(\textbf{Figure 1a, b} ). The final assembly of \emph{S.
furcifera}(`Surf\_2022') comprised 642,018,011 base pairs, with a
scaffold N50 of 41,905,432 base pairs and N90 of 21,229,120 base pairs.
`Surf\_2022' consisted of 15 chromosome-level scaffolds (with lengths
ranging from 13.80 to 74.02 Mb, total 621,787,945 base pairs) covering
95.04\% of the total sequences, together with an additional 624 small
scaffolds (20,230,066 base pairs) making up the remaining 4.96\% of the
genome (\textbf{Table 1} ). The GC content of `Surf\_2022' was 34.25\%.

To date, three versions of genome assemblies have been completed,
including the first version `Surf\_2017' based on the Illumina sequence
from our research group and two other versions of the chromosomal-level
assemblies `OGS' and `ASM1714138v1'. In terms of genome size, the
currently assembled reference genomes `Surf\_2022' (642 Mb) and `OGS'
(656 Mb) were close to the results of genome size assessed by BUSCO,
while the `ASM1714138v1' (561 Mb) genome differed by nearly 90 Mb. There
were also partial gap regions in the genome, with `Surf\_2022'
containing 588 gaps, which was lower than the 777 gaps in `ASM1714138v1'
and much lower than the 26609 gaps in `OGS' (\textbf{Table S3} ). These
results indicated that the `Surf\_2022' genome had higher integrity than
all the previous versions. Genome completeness was evaluated using BUSCO
assessment and the matching rate with the SRR3954848 genomic library.
The BUSCO for `Surf\_2022' reached 98.2\%, with 95.1\% single-copy
genes, 3.1\% double-copy, 1\% fragmentation, and 0.8\% deletion genes
(\textbf{Figure 1c} ). The BUSCO completeness of `ASM1714138v1' and
`OGS' was 95.9\% and 95.2\%, respectively, which was lower than that of
`Surf\_2022' (98.20\%). The above statistics and evaluation of genome
quality verified the highest completeness and accuracy of our genome
assembly at the chromosome level among all four versions of the genome.

\textbf{Genome annotation}

The genome landscape was generated to show the chromosome-level genome
characteristics of \emph{S. furcifera} (\textbf{Figure 1d} ). From the
assembled genome of \emph{S. furcifera} , a total of 24,807 genes were
predicted, of which 24,669 protein-coding genes were acquired
(\textbf{Table 1} ). Gene function annotation was assigned on the basis
of sequence and domain conservation. A total of 14535 genes shared
homology with known genes and were functionally annotated through public
databases. The gene number of annotations was 12,167 by GO, 11,406 by
Pfam, and 5,064 by KEGG. In the \emph{S. furcifera} genome, 401 KEGG
pathways were identified, and several genes in organismal systems,
metabolism, genetic information processing, environmental information
processing and cellular pathways were significantly enriched, such as
endocrine system and immune system in organismal systems pathways and
signaling transduction in environmental information processing pathways
(\textbf{Figure S2a} ). In addition, 19,629 GO terms were identified,
including 10,164 genes involved in 14,311 biological processes, 10,768
genes possessing 3,540 kinds of molecular functions, and 8,991 genes
containing 1,778 cellular components that were significantly enriched,
such as binding and protein binding in the molecular function category
(\textbf{Figure S2b} ).

Additionally, we identified 1,211 lncRNAs across the whole genome, the
majority (74.5\%) of which were longer than 1,000 bp (\textbf{Table S4}
). The developmental expression profiles of lncRNAs at eight different
developmental stages are shown in heatmaps (\textbf{Figure S3a} ). It is
worth noting that the embryo stage had more lncRNAs showing upregulated
expression than the other stages, suggesting that the effect of lncRNAs
may be more pronounced at these three developmental stages. Furthermore,
the different lncRNAs revealed expression differences at other
developmental stages, especially in 1\textsuperscript{st} and
5\textsuperscript{th} instars, 5-day-old adults and 10-day-old adults.
Furthermore, 7,595 circRNAs were also identified, and the different
circRNAs revealed significant expression differences at all eight
developmental stages (\textbf{Figure S3b} ), suggesting that circRNAs
may have important effects on WBPH development.

Finally, the repeat sequences were also estimated to account for 46.40\%
(297.88 Mb) of the \emph{S. furcifera} assembly, among which 24.01\%
were uncharacterized repeats. In the remaining parts, DNA and long
interspersed nuclear elements (LINEs) were predominant (40.33 Mb, 6.28\%
and 41.05 Mb, 6.39\%, respectively), with long terminal repeats (LTRs)
being the most abundant subtype, representing 4.54\% (29.14 Mb) of the
genome assembly. The results of repeat annotations are summarized
in\textbf{Table 2} .

\textbf{Comparative genomics}

Identification of orthologous genes and phylogenetic analysis

Totally, 860 strict single-copy orthologs were used to perform
phylogenetic analysis between \emph{S. furcifera} and 16 other species
(\textbf{Figure 2a, b} ). As expected, the phylogenetic analysis
identified the closest relationship of \emph{S. furcifera} to \emph{L.
striatellus} , which further formed a sister taxon with \emph{N. lugens}
, supporting the closer relative relationships in the same family of
Hemiptera compared with other insect species. Divergence time analysis
showed that Hemiptera (hemimetabolous insects) diverged from
Hymenoptera, Coleoptera, Lepidoptera and Diptera (holometabolous
insects) approximately 392 million years ago (Mya),\emph{S. furcifera}
diverged with\emph{L. striatellus} from the most recent common ancestor
approximately 46 Mya, and \emph{N. lugens} diverged from \emph{L.
striatellus} and\emph{S. furcifera} approximately 74 Mya (\textbf{Figure
2a} ). In addition, a total of 244,544 genes were found to have
orthologous genes shared between at least two arthropod species. The
largest number of genes of\emph{S. furcifera} (19,799) was found to be
orthologous to 22,029 genes of \emph{L. striatellus} , followed by
17,946 genes of \emph{S. furcifera,} which were orthologous to 19,818
genes of \emph{N. lugens}among the studied insects, further supporting
the close relationships in the same insect order. In addition, 2,516
genes were identified to be specific to the \emph{S. furcifera} genome
when compared with the other 16 genomes (\textbf{Figure 2a} ).

Gene family expansion and contraction analysis

The expansion and contraction of gene families may indicate functional
adaptation and evolution. According to the CAF\selectlanguage{ngerman}É analysis, we found that
1,190 and 1,662 gene families were expanded and contracted in \emph{S.
furcifera} , respectively (\textbf{Figure 2a} ). Moreover, a total of
173 rapidly evolving gene families (\emph{P} value \textless{} 0.05)
(148 expanded and 25 contracted families) were found in \emph{S.
furcifera} . Functional enrichment showed that these rapidly expanded
gene families were associated with multiple immune-related pathways,
such as the NF-kappa B signaling pathway, ubiquitin-mediated
proteolysis, nucleocytoplasmic transport, and metabolism of xenobiotics
by cytochrome P450 (\textbf{Figure 2b} ). Additionally, we found that
the three rice planthoppers we most concerned in Delphacidae of
Hemiptera had undergone rapid co-expansions and co-contractions of
conserved genes and gene families during evolution through comparative
analyses (\textbf{Figure 2c} ). In total, 79 gene families revealed
rapid co-expansion, including 1246 genes of brown planthopper (BPH), 634
genes of small brown planthopper (SBPH), and 1741 genes of white-backed
planthopper (WBPH). In particular, as BPH and WBPH possess long-distance
migration ability, while SBPH does not, we also paid attention to the
gene families that were co-expanded in BPH and WBPH but contracted in
SBPH. We found that three orthogroups (OG0000044, OG0000047, OG0000063)
were rapidly co-expanded in BPH and WBPH but rapidly contracted in SBPH,
including 87 genes in WBPH, 91 genes in BPH and 20 genes in SBPH, of
which most genes in OG0000047 shared a common domain, PF13843, which was
annotated as piggyBac transposable element-derived protein 4 (PGBD4).

Furthermore, we performed positive selection analysis on the genes
associated with immunity among three rice planthoppers in Delphacidae of
Hemiptera. A total of 77 genes associated with immunity were identified
to undergo selection pressure among the three rice planthoppers and are
thus designated positively selected genes (PSGs) in \textbf{Table S5} .
We found that several PSGs are core components in classic immune
pathways, such as the \emph{bsk} , \emph{kay} , and \emph{Uev1A} genes
in the Imd signaling pathway; the \emph{pll,} \emph{spz} , \emph{Toll-7}
and\emph{Toll-9} genes in the Toll-NF-\selectlanguage{greek}κ\selectlanguage{english}B signaling pathway; and the
GNBP1 and GNBP3 genes in the extracellular sp\selectlanguage{ngerman}ätzle activating pathway.

Immune gene identification and RNAi effects of the genes in classic
immune pathways

As described above, the genes in several immune-related pathways had
undergone rapid evolution and positive selection; therefore, we further
performed genome-wide identification and annotation of the immune genes
in \emph{S. furcifera} . In total, 749 genes in classic immune pathways
were identified and annotated in the \emph{S. furcifera} genome,
including 84 genes in the RNAi pathway, 98 genes in the Toll pathway,
162 genes in the IMD pathway, 51 genes in the Jak-STAT pathway, 265
genes in the autophagy pathway, 72 genes in the phagocytosis pathway,
and 17 genes in the apoptosis pathway (\textbf{Table S6} ). In addition,
the localizations of these genes were also visualized on the respective
chromosomes (\textbf{Figure S4} ). As multiple immune-related genes in
NF-kappa B signaling pathways were detected to undergo rapid evolution
and positive selection, a simplified diagram was generated to provide an
overview of the core components (\textbf{Figure 3a, b} ). Notably, we
found that the FADD (Fas-associated death domain) gene in the
well-conserved Imd immune pathway is missing in all species of Hemiptera
except for \emph{Cimex lectularius} compared with all nine other insect
species in Diptera, Lepidoptera, Coleoptera, Hymenoptera and
Phthiraptera that were included in this study (\textbf{Figure 3b, Figure
S5} ).

Thereafter, we chose 14 representative immune genes that play important
roles in classic immune pathways for antiviral functional validation.
First, the viruliferous rates obtained after silencing of the genes were
compared (\textbf{Table S7} ). In the control dsGFP group for the immune
pathway, 36.36\% of the treated nymphs were viruliferous at 24 hours
after dsRNA injection. Knockdown of all target genes in immune pathways
caused an increase in the viruliferous rate (41.18\% to 70.37\%). Then,
qPCR analysis revealed that each dsRNA efficiently suppressed the
transcript levels of its target genes (\textbf{Figure 3} ). Knockdown of
the three genes (\emph{Toll} , \emph{Dorsal} and \emph{Dif} ) resulted
in an increase in the replication of the SRBSDV compared with the
control, especially the significant effects from the knockdown of the
\emph{Dorsal} and\emph{Dif} genes (\emph{P} value \textless{} 0.001 by
two-sided T test) (\textbf{Figure 3c} ). In addition, the knockdown of
the two genes (\emph{IMD} and \emph{Relish} ) in the IMD pathway also
led to some extent but not a significant increase in SRBSDV replication
(\emph{P}value = 0.0752 and 0.1134, respectively, by two-sided T test).
However, the knockdown of the two genes (\emph{Hop} and \emph{STAT} ) in
the Jak-STAT pathway showed nonsignificant effects on SRBSDV replication
(\emph{P} value = 0.5295 and 0.0763, respectively, by two-sided T test).
In the control dsGFP group for autophagy, phagocytosis, and apoptosis
pathways, 28\% of the treated nymphs were viruliferous at 24 hours after
dsRNA injection. Knockdown of most target genes in these pathways caused
a slight increase in the viruliferous rate, while knockdown
of\emph{Eater} in the phagocytosis pathway and \emph{Diap1} in the
apoptosis pathway led to an increase in the viruliferous rate to 41.67\%
and 57.14\%, respectively, and knockdown of \emph{Atg13} in the
autophagy pathway significantly decreased the viruliferous rate to 13\%
(\textbf{Table S7} ). Knockdown of \emph{Eater} in the phagocytosis
pathway and \emph{Diap1} in the apoptosis pathway resulted in an
increase but not significant replication of SRBSDV compared with the
control (\emph{P} value = 0.8110 and 0.6411, respectively, by two-sided
T test). However, knockdown of\emph{Atg13}and \emph{Vps15} in the
autophagy pathway and \emph{Dcp1} in the apoptosis pathway resulted in a
significant decrease in the replication of SRBSDV compared with the
control (\emph{P} value \textless{} 0.0001 by two-sided T test)
(\textbf{Figure 3d} ). Taken together, these results indicated that the
\emph{Dorsal} and \emph{Dif} genes in the Toll immune pathway worked
well with the antiviral effects, while the\emph{Atg13} and
\emph{Vps15}genes in the autophagy pathway and \emph{Dcp1} in the
apoptosis pathway may play a role in assisting of the virus replication
during the SRBSDV infection.

\textbf{Population genomics}

Phylogenetic analyses of the WBPH population

\emph{S. furcifera} is a migratory insect that can migrate from southern
Asia, such as Vietnam and the southern part of China, to the northern
part of China. For a better understanding of the population
differentiation of the WBPH, we selected 44 individuals to represent 12
geographic populations of the WBPH in China and Vietnam for genome
resequencing, including 9 accessions from Vietnam and 35 accessions from
different parts of China (\textbf{Figure S6} ). A total of 43.5 Mb of
SNPs was identified in 44 individuals, with 2.6 Mb of SNPs identified
from each WBPH sample. Only 1,059,267 SNPs with no missing data were
retained for further analysis, and the density of SNPs reached
approximately 2 SNPs per kb, with 606,343 (57\%) from gene regions
(1026491) and 97\% from noncoding regions.

Population structure and gene flow of the WBPH population

To investigate the genetic divergence among WBPH populations, we
performed phylogenetic analysis on 12 geographic populations using the
genome-wide 465,882 SNPs (to ensure that adjacent SNPs were not linked,
the squared correlation (r\textsuperscript{2}) greater than 0.2 was
filtered out, from 1,059,267 to 465,882 SNPs) (\textbf{Figure 4a} ).
Among the 12 populations, the geographically nearby populations were
more closely related, such as HNQH (Hainan Qionghai), HNLS (Hainan
Lingshui) and HNBT (Hainan Baoting) from Hainan, China. A minimum
spanning network was constructed to show the genetic relatedness among
the multillocus genotypes based on the dissimilarity values
(\textbf{Figure 4b} ). We found that the samples from HNBT (T2, T3) and
VNTP (D2, D3) were distantly related to all the other samples. Moreover,
sample F3 from Bac Giang, northern Vietnam, was clustered in the center
of the network and closely related to the other genotypes with one or
more steps of genetic variation, indicating that northern Vietnam may
contribute to the yearly immigrations of the WBPH. The TreeMix software
indicated similar relationships between populations to those revealed in
the phylogenetic tree (\textbf{Figure 4c} ), and gene flow from south to
north was indicated, such as VNTP to GDHF, YNCJ and YNMS, VNBJ to YNMS,
and HNQH to YNCJ. Similarly, we also found reverse directions of gene
flow, such as YNSD to HNQH and VNND, GXYZ to VNND, and GDGZ to VNND. We
found that samples northeastern Vietnam and its adjacent locations in
Yunnan and Guangdong of China have varying degrees of gene flow. This is
consistent with our analysis of minimum networks (\textbf{Figure 4b} ),
suggesting the possible migration of WBPH populations from Vietnam to
China. Similarly, there is a certain degree of gene flow from Yunnan to
Vietnam, indicating that bidirectional gene exchange between China and
Vietnam may be one of the factors causing the absence of population
structure among the 12 populations.

To identify the population structure, we used the non-model clustering
method DAPC to generate any possible population structures using the
genome-wide 465,882 SNPs in our data. First, DAPC was performed for the
44 accessions without predefining the populations, and plotting the
value of BIC against the number of clusters (\textbf{Figure S7a} )
indicated that a K of 2 with an elbow was the optimal cluster number.
Interestingly, when K was increased to 6, new divergence only appeared
in the samples from Hainan (\textbf{Figure S7b} ). Second, we also
performed DAPC for the 44 accessions based on the 12 predefined
populations based on their geographies. Two kinds of DAPC analysis both
confirmed the results from the phylogenetic analysis and supported
nonclear population structure among the other 12 geographic populations,
while some samples from local regions revealed a clear separation
(\textbf{Figure 4d} ). In addition, the genotype composition plot
representing the group assignment probability of individuals to several
groups revealed population admixtures among most geographic populations
except for some samples in VNTP and HNBT, indicating the possible local
population formed in some regions, which coincided with the results from
the phylogenetic analysis (\textbf{Figure 4e} ). Taken together, our
analysis of genome-wide SNPs collectively supports that there was no
significant population structure of 12 geographic populations, but some
samples in certain locations were clearly separated from others,
suggesting that the migration of WBPH removes the gene variation among
most of the geographic populations collected in this study, while
population divergence still exists in some local regions.

Selective sweep analysis

To identify genomic regions in the WBPH showing signatures of selection
during evolution and adaptation, selective sweeps were conducted using
the resequencing data. Since no significant population structure was
detected in the aforementioned analysis, we merged all 44 accessions and
assumed them to be random samples from a randomly mating population. The
CLR (composite likelihood ratio) scores along scaffolds were calculated
and plotted in \textbf{Figure 5a} . Scanning the genome in 10-kb
windows, the empirical cutoff value was set to 1\% for the outlier
regions with significant deviations from neutrality, and 1.19\%
(2926/24,669) of genes within these regions were identified to be under
selection (\textbf{Figure 5a} ). Among them, 2432 genes have homologs in
other species, while 309 and 185 are single-copy and multiple-copy
specific genes of the WBPH, respectively.

KEGG analyses demonstrated that the subset of genes with the strongest
selective sweep signal was enriched in the pathways of energy metabolism
(vitamin digestion and absorption, arachidonic acid metabolism, linoleic
acid metabolism) and immunity (viral life cycle, phagosome, NF-kappa B
signaling pathway) (\textbf{Figure 5b} ). Gene Ontology (GO) analysis
identified a significant overrepresentation of genes involved in
biological processes that contribute to morphogenesis and metamorphosis
during insect development, including appendage morphogenesis,
metamorphosis, postembryonic animal morphogenesis and instar larval or
pupal morphogenesis. These processes were related to energy metabolism,
immunity and morphogenesis (\textbf{Figure 5b, c} ), which are very
important for WBPH environmental adaptation and long-distance migration.

Furthermore, we investigated the genes with the top strong selective
sweep signals across the genome (\textbf{Figure 5a} ), which were mainly
associated with morphogenesis, metamorphosis and other crucial processes
during insect development (\emph{Src64B,Galphao} , \emph{gbb} ,
\emph{ena} , \emph{kis} , \emph{mirr} ,\emph{emc} , \emph{Pms2} ,
\emph{mtt} , \emph{fuss} , \emph{Trf2} ,\emph{dac} , \emph{AnxB9} ,
\emph{Dcr-1} , \emph{wg, etc.} ), and immunity (\emph{psidin} in the
phagocytosis pathway, \emph{Toll-7} in the Toll pathway, \emph{IKKbeta}
in the Imd pathway, \emph{pkc53E} in the Hh (Hedgehog) signaling
pathway, which regulates the Duox-dependent gut immune response to
bacterial uracil,\emph{Atg1} and \emph{ras85D} in the autophagy
pathway). The genes \emph{fuss} , \emph{Trf2} , \emph{dac} ,\emph{AnxB9}
, \emph{Dcr-1} , and \emph{wg} play crucial roles in wing morphogenesis,
of which \emph{wg} is the core member in the wingless signaling pathway
that has been identified as an important gene controlling wing
development. Notably, \emph{IKKbeta} , \emph{Toll-7} and\emph{Atg1} are
core components in their respective immune pathways. Taken together,
selective sweep analysis provided new insights into the crucial roles of
the related genes involved in the morphogenesis and immunity of \emph{S.
furcifera} .

\textbf{DISCUSSION}

White-backed planthopper has been a huge threat to rice crops in China
and Southeast Asian countries (Liu, 2010), and it causes a double blow
for rice production by transmitting SRBSDV in a persistent and
propagative manner. Genome resources can provide important insights into
the complicated relationships between the WBPH and SRBSDV. To date,
three versions of genome assemblies of \emph{S. furcifera} have been
reported in public (Ma et al. 2021; Wang et al., 2017). However, all
three versions have their drawbacks, either low quality or insufficient
completeness. In the present study, a chromosome-level genome assembly
and annotation was completed by combining genomic reads from Nanopore
sequencing and transcriptomic reads from PacBio sequencing. The
chromosomal-level genome assembly not only reveals high accuracy but
also more insights into chromosome evolution, providing more genetic
resources to discover key genes or targets.

The \emph{S. furcifera} genome has a G+C content of 34.25\%, which is
slightly higher than that of our first version of the genome assembly
(31.6\%); compared to other hemipteran species, this is similar to that
of the small brown planthopper \emph{L. striatellus} (34.54\%) (Zhu et
al., 2017) and brown planthopper \emph{N. lugens} (34.6\%) (Xue et al.,
2014). For genome annotation, in the first version, a total of 21,254
protein-coding genes were annotated (Wang et al., 2017). In the third
version, a total of 15,929 protein-coding genes were annotated (Ma et
al., 2021). In the present study, 24,669 protein-coding genes were
annotated. In addition, ncRNAs constitute a major component of the
transcriptome, playing crucial roles in the transcriptional and
posttranscriptional regulation of gene expression (Mattick \& Makunin,
2006). In fruit flies, we have identified one lncRNA, VINR, that
activates noncanonical innate immune signaling upon detection of the
dsRNA-binding viral suppressors of RNAi (VSR) of Drosophila C virus
(DCV) by our research group (Zhang et al., 2020). For WBPH, ncRNAs have
been genome-widely characterized from the first version of the WBPH
genome based on next-generation sequencing, including 171 microRNAs
(Chang et al., 2018; Chang et al., 2016) and 1861 lncRNAs (Chang, Ajayi,
Guo, \& Wu, 2020). However, functional validation of these ncRNAs is
still deficiency in WBPH. Therefore, the 1,211 lncRNAs and 7595 circRNAs
identified in the current WBPH genome plus the sufficient experience and
technology from the study of fruit flies promote our depth study of the
lncRNAs and circRNAs functioning during antiviral and developmental
processes of WBPH. Additionally, we obtained high-quality transcriptomic
data from SRBSDV-infected and SRBSDV-free WBPHs and found some potential
lncRNAs that may play important roles during antiviral processes, which
will be our next step to deeply detect the immune mechanisms with
large-scale experimental validation.

Comparative genomic analysis was performed to obtain an overview of the
genetic evolution among different insect species. Notably, we found that
one orthogroup (OG0000047) shared a common domain, PF13843 (PGBD4),
which was rapidly co-expanded in BPH and WBPH (long-distance migratory
insects) but rapidly contracted in SBPH (non-long-distance migratory
insect in most situations). PGBD4 has been reported as a target of a
lncRNA to decrease the expression of canonical wing development-related
genes and ultimately regulate wing development in the brown citrus
aphid\emph{Aphis citricidus} (Shang et al., 2021). The development of
wings is crucial for insects to adapt to environmental stresses,
allowing winged morphs to emigrate from unfavorable environments to new
habitats (Zhang et al., 2019). This phenomenon promotes our future study
of any possible function of the PGBD4 gene in wing development,
especially for the migratory insects WBPH and BPH. Moreover, 148 rapidly
expanded gene families in \emph{S. furcifera} were identified that were
associated with multiple immune-related pathways (\textbf{Figure 2b} ).
Furthermore, we found that 77 genes associated with immunity underwent
selection pressure among the three rice planthoppers (\textbf{Table S5}
), of which several genes are core components in classic immune
pathways. Together, these results indicate that \emph{S. furcifera}
might have enhanced adaptation to environmental stress. All these
findings promote our further study of the immune genes of \emph{S.
furcifera} . Notably, we found that the core component FADD gene in the
well-conserved Imd immune pathway is missing in most of the Hemiptera
species used in this study (\textbf{Figure 3b} ), which is deserved to
determine whether there were substituted genes in this pathway or any
novel mechanism working here. Fourteen representative immune genes were
selected for functional validation, demonstrating the antiviral effects
of \emph{Dorsal} and\emph{Dif} genes and assisting virus replication of
\emph{Atg13} ,\emph{Vps15} and \emph{Dcp1} genes during the SRBSDV
infection for WBPH.\emph{Dorsal} and \emph{Dif} are two very similar
NF-kB homologs, of which \emph{Dorsal} is an invertebrate NF-\selectlanguage{greek}κ\selectlanguage{english}B
transcription factor that is essential for the regulation of immune
responses, and \emph{Dif} is a\emph{Dorsal} -related immunity factor (Ip
et al., 1993; Manfruelli, Reichhart, Steward, Hoffmann, \& Lemaitre,
1999). \emph{Dorsal} has been demonstrated to be required for \emph{D.
melanogaster} to resist oral infection with Drosophila C virus (Ferreira
\selectlanguage{ngerman}Á et al., 2014). In \emph{S. furcifera} , 348 immune-related genes were
identified from homology-based analysis, and the RNAi pathway was
suggested as the major antiviral response based on differential
expression analysis (Wang et al., 2016). However, these conclusions
remain only in terms of data analysis based on the next generation
sequencing data with a lower completeness and accuracy, and far less is
known about immune genes in other immune pathways. Therefore, the more
accurate and systematic identification of immune genes based on
high-quality genomic and transcriptomic data from third-generation
sequencing and RNAi experiments of silencing the representative immune
genes in this study provided novel insights for understanding the immune
mechanisms. Subsequent experiments are needed for the deep detection of
the immune interaction of the key genes in the antiviral process.

Genome-wide SNPs were identified by resequencing 44 individuals from 12
geographic populations collected in China and Vietnam. The DAPC and
phylogenetic analysis using the high density of SNP information revealed
nonclear genetic structure of the 12 geographic populations, while some
local populations were indicated in Hainan, China. The genetic structure
of WBPH has been investigated previously by several studies using
microsatellites or mitochondrial gene markers (Liu, Gui, Li, 2013;
Matsumoto et al., 2013; Mun et al., 1999; Sun, Jiang, Wang, Hong, 2014;
Xie, Guo, Jin, Wang, 2014). However, there were no reports showing that
the distinct WBPH population in Hainan, either because of the deficiency
of sample collections in the local regions from Hainan or less robust
markers used in these studies. Our study using high-density SNPs clearly
demonstrates that significant differentiation existed in three samples
from the HNBT population. Hainan Island is one of the largest islands
with a complex landscape in China and has three main mountains, i.e.,
Wuzhishan (1867 m a.s.l.), Yinggeling (1811 m a.s.l.), and Mihouling
(1654 m a.s.l.) peak at the center, which divides the island into a
north plain and a south hilly land (Liu et al., 2022). These mountains
provide a possible barrier to limit migrations of the WBPH heading
north, because the weak gene exchange with other locations may shape a
single genetic structure in HNBT, which has some parts separated by
Wuzhishan from other locations in northern Hainan, such as HNQH.
Additionally, other geographic populations revealed no clear genetic
differentiation, which is also consistent with other pest populations
separated by large geographic distances, such as \emph{N. lugens}
(Matsumoto et al., 2013), \emph{L. striatellus} (Han et al., 2020) and
\emph{Plutella xylostella} (Endersby et al., 2006). The weak
phylogeographic structure between these Asian populations from
relatively remote geographic origins was probably caused by the high
dispersal ability of the WBPH that allowed genetic mixing and may also
be attributed to the East Asian subtropical monsoon facilitating the
WBPH migration and promoting the gene flow between these Asian
populations (Kwon et al., 2018). As the gene flow analysis shows in
\textbf{Figure 4c} , the VNTB and VNBG in northern Vietnam and HNQH in
China had high migration rates to Guangdong and Yunnan in China, which
is consistent with previous reports that both the central and northern
parts of Vietnam and southern Hainan Province in China are source areas
for early immigration to Guangdong and Guangxi (Zhai et al., 2011).
Meanwhile, bidirectional gene flow was also revealed in the tree,
indicating the possible influence of different monsoons. However, to
obtain clearer information on population evolution and migration, we
need more sample collections in a wider geographic range in the future.

Sweep selective analysis was also performed, and the top strong
selective sweep signals across the genome (\textbf{Figure 5a} ) were
mainly associated with morphogenesis, metamorphosis and other crucial
processes during insect development (\emph{Src64B,} \emph{Galphao}
,\emph{gbb} , \emph{ena} , \emph{kis} , \emph{mirr} , \emph{emc}
,\emph{Pms2} , \emph{mtt} , \emph{fuss} , \emph{Trf2} , \emph{dac}
,\emph{AnxB9} , \emph{Dcr-1} , \emph{wg, etc.} ), and immunity
(\emph{psidin} , \emph{Toll-7} , \emph{IKKbeta} , \emph{pkc53E}
,\emph{Atg1} and \emph{ras85D} ). Notably, \emph{IKKbeta} ,
\emph{Toll-7}and \emph{Atg1} are core components in their respective
immune pathways. Combined with the results from the comparative genomic
analysis, we gained a consistent conclusion that the immune-related
genes of \emph{S. furcifera} have undergone rapid evolution and should
be given more attention in future experimental research. Moreover, the
genes\emph{fuss} , \emph{Trf2} , \emph{dac} , \emph{AnxB9} ,
\emph{Dcr-1} , and\emph{wg} play crucial roles in wing morphogenesis
(Akai, Ohsawa, Sando, \& Igaki, 2021; Fischer et al., 2012; Klein, 2001;
O'Keefe et al., 2012), of which \emph{wg} is the core member in the
wingless signaling pathway that has been identified as an important gene
controlling wing development in the fruit fly \emph{D. melanogaster}
(Klein, 2001). In our research group, 45 genes in the module directly
interacting with the IIS-PI3K-AktFOXO pathway have been regarded as
potential insulin pathway-mediated wing dimorphism-related genes, 5 of
which have been experimentally verified for insect wing size regulation
(Gao et al., 2019). Therefore, these genes that have undergone selection
pressure and associated with wing development will also have
experimental verification in the future and will provide more
information for the study of the evolutionary adaptation of the
planthopper.

The chromosome-level genome assembly of\emph{S. furcifera} would provide
a valuable resource for discovering the key protein-coding genes and
noncoding RNAs, identifying the genetic evolution and population
structure, further detecting the interaction mechanisms behind SRBSDV
and WBPH, and developing an effective pest management system.

\textbf{Table 1} Summary of genome sequencing, assembly, and annotation.\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
\textbf{Genome Assembly and Annotations} &
\textbf{Statistics}\tabularnewline
\midrule
\endhead
\textbf{Assembly} &\tabularnewline
Number of contigs & 639\tabularnewline
Contig N50 (bp) & 41,905,432\tabularnewline
Contig N75 (bp) & 1,898,492\tabularnewline
Maximum contig size (bp) & 19,926,834\tabularnewline
Genome size (bp) & 642,018,011\tabularnewline
Number of chromosomes & 15\tabularnewline
Number of scaffolds & 624\tabularnewline
Length of chromosomes (bp) &
13,806,739\textasciitilde{}74,022,170\tabularnewline
GC content (\%) & 34.25\tabularnewline
\textbf{Annotations} &\tabularnewline
Total number of genes & 24,807\tabularnewline
Number of GO annotation & 12,167\tabularnewline
Number of Pfam annotation & 11,406\tabularnewline
Number of KEGG annotation & 5,064\tabularnewline
Number of protein-coding genes & 24,669\tabularnewline
Number of genes identified in classic immune signaling pathways &
749\tabularnewline
Number of lncRNAs & 1121\tabularnewline
Number of circRNAs & 7595\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2} Repeats in the\emph{Sogatella furcifera} genome
assembly.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
Class & Number & Length (bp) & Ratio of genome (\%)\tabularnewline
\midrule
\endhead
DNA & 149371 & 40327827 & 6.281416769\tabularnewline
LINEs & 119061 & 41046187 & 6.39330771\tabularnewline
SINEs & 35546 & 5068824 & 0.789514299\tabularnewline
LTRs & 73330 & 29140429 & 4.538880296\tabularnewline
Others & 335098 & 28122336 & 4.380303281\tabularnewline
Unknown & 718792 & 154170486 & 24.01342071\tabularnewline
Total & 1431198 & 297876089 & 46.39684306\tabularnewline
\bottomrule
\end{longtable}

Abbreviations: DNA, DNA elements; LINEs, long interspersed nuclear
elements; SINEs, short interspersed nuclear elements; LTRs, long
terminal repeats; Others, Small RNA, satellites and simple repeats.

\textbf{ACKNOWLEDGEMENTS}

We are grateful to Professor Chuanxi Zhang for supplying the genome
sequence of \emph{Nilaparvata lugens} . This work was supported by the
National Natural Science Foundation of China (grant 32030087). The
funders played no role in the study design, data collection and
analysis, decision to publish, or preparation of the manuscript.

\textbf{DATA AVAILABILITY STATEMENT}

The Whole genome project of WBPH has been deposited at the public NCBI
Nucleotide Database under BioProject (PRJNA923506).

\textbf{AUTHOR CONTRIBUTIONS}

Ting Cui: Formal analysis, Software, Investigation, Writing-review and
editing; Qing Bai: Formal analysis, Writing-original draft;
Writing-review and editing; Weiqiang Yu: Formal analysis, Software,
Methodology; Dongyang Guo: Methodology, Investigation; Yawen Ban:
Methodology; Kun Chen: Methodology; Ali Raza, Writing-review and
editing; Guohui Zhou: Sample collection; Qingfa Wu: Project
administration, Supervision, Sample Collection, Funding acquisition,
Writing-review and editing.

\textbf{ORCID}

\emph{Qingfa Wu} https://orcid.org/0000-0001-5420-6213

\emph{Qing Bai} https://orcid.org/0000-0002-3505-4426

\textbf{REFERENCES}

Akai, N., Ohsawa, S., Sando, Y., \& Igaki, T. (2021). Epithelial
cell-turnover ensures robust coordination of tissue growth in Drosophila
ribosomal protein mutants. \emph{PLoS genetics} ,~\emph{17} (1),
e1009300. https://doi.org/10.1371/journal.pgen.1009300

Bashir, A., Klammer, A. A., Robins, W. P., Chin, C.-S., Webster, D.,
Paxinos, E., \ldots{} \& Schadt, E. E. (2012). A hybrid approach for the
automated finishing of bacterial genomes. \emph{Nature Biotechnology,
30} (7), 701--707. https://doi.org/10.1038/nbt.2288

Burton, J. N., Adey, A., Patwardhan, R. P., Qiu, R., Kitzman, J. O., \&
Shendure, J. (2013). Chromosome-scale scaffolding of de novo genome
assemblies based on chromatin interactions. \emph{Nature Biotechnology,
31} (12), 1119--1125. https://doi.org/10.1038/nbt.2727

Cantarel, B. L., Korf, I., Robb, S. M., Parra, G., Ross, E., Moore, B.,
\ldots{} \& Yandell, M. (2008). MAKER: an easy-to-use annotation
pipeline designed for emerging model organism genomes. \emph{Genome
Research, 18} (1), 188--196. https://doi.org/10.1101/gr.6743907

Chang, Z.-X., Ajayi, O. E., Guo, D.-Y., \& Wu, Q.-F. (2020). Genome-wide
characterization and developmental expression profiling of long
noncoding RNAs in \emph{Sogatella furcifera} . \emph{Insect Science,27}
(5), 987--997. https://doi.org/https://doi.org/10.1111/1744-7917.12707

Chang, Z.-X., Akinyemi, I. A., Guo, D.-Y., \& Wu, Q.-F. (2018).
Characterization and comparative analysis of microRNAs in the rice
pest\emph{Sogatella furcifera} . \emph{PLoS One, 13} (9), e0204517.
https://doi.org/10.1371/journal.pone.0204517

Chang, Z.-X., Tang, N., Wang, L., Zhang, L.-Q., Akinyemi, I. A., \& Wu,
Q.-F. (2016). Identification and characterization of microRNAs in the
white-backed planthopper, \emph{Sogatella furcifera} . \emph{Insect
Science, 23} (3), 452--468. https://doi.org/10.1111/1744-7917.12343

Cook, D. E., Valle-Inclan, J. E., Pajoro, A., Rovenich, H., Thomma, B.
P. H. J., \& Faino, L. (2018). Long-Read Annotation: Automated
Eukaryotic Genome annotation based on long-read cDNA
sequencing.\emph{Plant Physiology, 179} (1), 38--54.
https://doi.org/10.1104/pp.18.00848

Darriba, D., Taboada, G. L., Doallo, R., \& Posada, D. (2011). ProtTest
3: fast selection of best-fit models of protein
evolution.\emph{Bioinformatics} , \emph{27} , 1164--1165.
https://doi.org/10.1093/bioinformatics/btr088

O'Keefe, D. D., Thomas, S. R., Bolin, K., Griggs, E., Edgar, B. A., \&
Buttitta, L. A. (2012). Combinatorial control of temporal gene
expression in the Drosophila wing by enhancers and core
promoters.\emph{BMC genomics} , \emph{13} (1), 1--15.

De Bie, T, Cristianini N, Demuth J. P., \& Hahn, M. W. (2006). CAFE: a
computational tool for the study of gene family
evolution.\emph{Bioinformatics} , \emph{22} , 1269--71.
https://doi.org/10.1093/bioinformatics/btl097

Dudchenko, O., Batra, S. S., Omer, A. D., Nyquist, S. K., Hoeger, M.,
Durand, N. C., \ldots{} \& Aiden, E. L. (2017). De novo assembly of
the\emph{Aedes aegypti} genome using Hi-C yields chromosome-length
scaffolds. \emph{Science, 356} (6333), 92--95.
https://doi.org/10.1126/science.aal3327

Dudchenko, O., Shamim, M. S., Batra, S. S., Durand, N. C., Musial, N.
T., Mostofa, R., \ldots{} \& Stamenova, E. J. B. (2018). The Juicebox
Assembly Tools module facilitates de novo assembly of mammalian genomes
with chromosome-length scaffolds for under \$1000. 254797.
https://doi.org/10.1101/254797

Durand, N. C., Shamim, M. S., Machol, I., Rao, S. S. P., Huntley, M. H.,
Lander, E. S., \& Aiden, E. L. (2016). Juicer provides a one-click
system for analyzing loop-resolution Hi-C experiments. \emph{Cell
Systems, 3} (1), 95--98. https://doi.org/10.1016/j.cels.2016.07.002

Emms, D. M., \& Kelly, S. (2019). OrthoFinder: phylogenetic orthology
inference for comparative genomics.\emph{Genome Biology, 20} (1), 238.
https://doi.org/10.1186/s13059-019-1832-y

Ferreira \selectlanguage{ngerman}Á, G., Naylor, H., Esteves, S. S., Pais, I. S., Martins, N. E.,
\& Teixeira, L. (2014). The Toll-dorsal pathway is required for
resistance to viral oral infection in \emph{Drosophila} . \emph{PLoS
Pathogen, 10} (12), e1004507.
https://doi.org/10.1371/journal.ppat.1004507

Fischer, S., Bayersdorfer, F., Harant, E., Reng, R., Arndt, S.,
Bosserhoff, A. K., \& Schneuwly, S. (2012). fussel (fuss)-A Negative
Regulator of BMP Signaling in \emph{Drosophila melanogaster} .
\emph{PLoS ONE} , \emph{7} (8): e42349.
https://doi.org/10.1371/journal.pone.0042349

Flynn, J. M., Hubley, R., Goubert, C., Rosen, J., Clark, A. G.,
Feschotte, C., Smit, A. F. (2020). RepeatModeler2 for automated genomic
discovery of transposable element families. \emph{Proceedings of the
National Academy of Sciences of the United States of America 117} (17),
9451--9457. https://doi.org/10.1073/pnas.1921046117

Gao, J., Zhao, B.-R., Zhang, H., You, Y.-L., Li, F., \& Wang, X.-W.
(2021). Interferon functional analog activates antiviral Jak/Stat
signaling through integrin in an arthropod. \emph{Cell Reports, 36}
(13), 109761. https://doi.org/10.1016/j.celrep.2021.109761

Gao, X.-L., Fu, Y.-T., Ajayi, O. E., Guo, D.-Y., Zhang, L.-Q., \& Wu,
Q.-F. (2019). Identification of genes underlying phenotypic plasticity
of wing size via insulin signaling pathway by network-based analysis
in\emph{Sogatella furcifera} . \emph{BMC Genomics, 20} (1), 396.
https://doi.org/10.1186/s12864-019-5793-z

Han, L., Zhang, J.-T., Wang, M.-M., Zhu, K.-X., \& Wang, X.-Y. (2020).
Mitochondrial DNA diversity and population structure of\emph{Laodelphax
striatellus}across a broad geographic area in China. \emph{Mitochondrial
DNA Part A, 31} (8), 346--354.
https://doi.org/10.1080/24701394.2020.1830075

Hoang, A. T., Zhang, H.-M., Yang, J., Chen, J.-P., Hébrard, E., Zhou,
G.-H., \ldots{} \& Cheng, J.-A. (2011). Identification,
characterization, and distribution of southern rice black-streaked dwarf
virus in Vietnam.\emph{Plant Disease} , \emph{95} (9), 1063--1069.
https://doi.org/10.1094/PDIS-07-10-0535

Hoff, K. J., Lange, S., Lomsadze, A., Borodovsky, M., \& Stanke, M.
(2015). BRAKER1: Unsupervised RNA-Seq-Based Genome Annotation with
GeneMark-ET and AUGUSTUS. \emph{Bioinformatics, 32} (5), 767--769.
https://doi.org/10.1093/bioinformatics/btv661

Hu, J., Fan, J.-P., Sun, Z.-Y., \& Liu, S.-L. (2019). NextPolish: a fast
and efficient genome polishing tool for long-read
assembly.\emph{Bioinformatics} , \emph{36} (7), 2253--2255.
https://doi.org/10.1093/bioinformatics/btz891

Huddleston, J., Ranade, S., Malig, M., Antonacci, F., Chaisson, M., Hon,
L., \ldots{} \& Eichler, E. E. (2014). Reconstructing complex regions of
genomes using long-read sequencing technology. \emph{Genome Research,
24} (4), 688--696. https://doi.org/10.1101/gr.168450.113

Ip, Y. T., Reach, M., Engstrom, Y., Kadalayil, L., Cai, H.,
González-Crespo, S., \ldots{} \& Levine, M. (1993). Dif, a
dorsal-related gene that mediates an immune response in Drosophila.
\emph{Cell} ,\emph{75} (4), 753--763.
https://doi.org/10.1016/0092-8674(93)90495-C

Jombart, T. (2008). adegenet: a R package for the multivariate analysis
of genetic markers. \emph{Bioinformatics} , \emph{24} , 1403--1405.
https://doi.org/10.1093/bioinformatics/btn129

Jombart, T., Devillard, S., \& Balloux, F. (2010). Discriminant
analysisi of principal componets: a new method for the analysis of
genetically structured populations. \emph{BMC Genetics} , \emph{11} ,
94. https://doi.org/10.1186/1471-2156-11-94

Jensen, L. J., Julien, P., Kuhn, M., von Mering, C., Muller, J., Doerks,
T., \& Bork, P. (2008). eggNOG: automated construction and annotation of
orthologous groups of genes. \emph{Nucleic Acids Research, 36} (Database
issue), D250--254. https://doi.org/10.1093/nar/gkm796

Jones, P., Binns, D., Chang, H.-Y., Fraser, M., Li, W., McAnulla, C.,
\ldots{} \& Hunter, S. (2014). InterProScan 5: genome-scale protein
function classification. \emph{Bioinformatics} , \emph{30} (9),
1236--1240. https://doi.org/10.1093/bioinformatics/btu031

Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M., \& Tanabe, M.
(2016). KEGG as a reference resource for gene and protein
annotation.\emph{Nucleic Acids Research} , \emph{44} , D457--462.
https://doi.org/10.1093/nar/gkv1070

Kingsolver, M. B., Huang, Z., \& Hardy, R. W. (2013). Insect antiviral
innate immunity: pathways, effectors, and connections.~\emph{Journal of
molecular biology} , \emph{425} , 4921--4936.
https://doi.org/10.1016/j.jmb.2013.10.006

Klein, T. (2001). Wing disc development in the fly: the early
stages.\emph{Current opinion in genetics \& development} , \emph{11}
(4), 470--475. https://doi.org/10.1016/S0959-437X(00)00219-7

Korf, I. (2004). Gene finding in novel genomes. \emph{BMC
bioinformatics} , \emph{5} (1), 1--9.
https://doi.org/10.1186/1471-2105-5-59

Kwon,~D. H.,~Jeong,~I. H.,~Lee,~S. W.,~\& Lee,~S. H.~2018.~Comparing the
genetic diversity and structure of indigenous Korean and Chinese
populations of~\emph{Laodelphax striatellus} ~fallén using mitochondrial
haplotypes. \emph{Journal of Asia-Pacific Entomology} , \emph{21} (1),
233--238. https://doi.org/10.1016/j.aspen.2017.12.004

Lin, H.-H., \& Liao, Y.-C. (2015). Evaluation and Validation of
Assembling Corrected PacBio Long Reads for Microbial Genome Completion
via Hybrid Approaches. \emph{PLoS One} , \emph{10} (12), e0144305.
https://doi.org/10.1371/journal.pone.0144305

Li, H. \& Durbin, R. (2009). Fast and accurate short read alignment with
Burrows\selectlanguage{english}-Wheeler transform. \emph{Bioinformatics} , \emph{25} ,
1754--1760. http://doi:.org/10.1093/bioinformatics/btp324.

Liu, J.-N., Gui, F.-R. \& Li, Z.-Y. (2010). Genetic diversity of the
planthopper, \emph{Sogatella furcifera} in the Greater Mekong Subregion
detected by intersimple sequence repeats (ISSR) markers. \emph{Journal
of Insect Science} , \emph{10} (1), 52.
https://doi.org/10.1673/031.010.5201

Liu, P., Zhang, X.-L, Huang, Q., Rizo, E. Z. C., Dumont, H. J., \& Han,
B.-P. (2022). Genetic diversity, lineage divergence, and demography
of\emph{Diaphanosoma dubium} (Crustacea: Cladocera) on Hainan Island,
China. \emph{Hydrobiologia} , \emph{849} (19), 4221--4239.
https://doi.org/10.1007/s10750-022-04975-w

Liu, W.-C., Liu, Y., Guo, R. (2010). Warning to outbreak of Southern
rice black-streaked dwarf virus. \emph{China Plant Protection}
,\emph{30} (3), 17--18.

Lomsadze, A., Burns, P. D., \& Borodovsky, M. (2014). Integration of
mapped RNA-Seq reads into automatic training of eukaryotic gene finding
algorithm. \emph{Nucleic Acids Research} , \emph{42} (15), e119.
https://doi.org/10.1093/nar/gku557

Ma, W.-H., Xu, L., Hua, H.-X., Chen, M.-Y., Guo, M.-J., He, K., \ldots{}
\& Li, F. (2021). Chromosomal-level genomes of three rice planthoppers
provide new insights into sex chromosome evolution. \emph{Molecular
Ecology Resources} , \emph{21} (1), 226--237.
https://doi.org/10.1111/1755-0998.13242

Manfruelli, P., Reichhart, J. M., Steward, R., Hoffmann, J. A., \&
Lemaitre, B. (1999). A mosaic analysis in \emph{Drosophila} fat body
cells of the control of antimicrobial peptide genes by the Rel proteins
Dorsal and DIF. \emph{The EMBO Journal} , \emph{18} (12), 3380--3391.
https://doi.org/10.1093/emboj/18.12.3380

Matsukura, K., Towata, T., Sakai, J., Onuki, M., Okuda, M., \&
Matsumura, M. (2013). Dynamics of Southern rice black-streaked dwarf
virus in rice and implication for virus
acquisition.\emph{Phytopathology} , \emph{103} (5), 509--512.
https://doi.org/10.1094/phyto-10-12-0261-r

Matsumoto, Y., Matsumura, M., Sanada-Morimura, S., Hirai, Y., Sato, Y.,
\& Noda, H. (2013). Mitochondrial cox sequences of \emph{Nilaparvata
lugens} and \emph{Sogatella furcifera} (Hemiptera, Delphacidae): low
specificity among Asian planthopper populations. \emph{Bulletin of
Entomological Research} , \emph{103} , 382--392.
https://doi.org/10.1017/S000748531200082X

Mattick, J. S., \& Makunin, I. V. (2006). Noncoding RNA. \emph{Human
molecular genetics} , \emph{15} (suppl\_1), R17--R29.
https://doi.org/10.1093/hmg/ddl046

McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K.,
Kernytsky, A., \ldots{} \& DePristo, M. A. (2010). The Genome Analysis
Toolkit: a MapReduce framework for analyzing next-generation DNA
sequencing data.~\emph{Genome research} ,~\emph{20} (9), 1297--1303.
https://doi.org/10.1101/gr.107524.110

Moriya, Y., Itoh, M., Okuda, S., Yoshizawa, A. C., \& Kanehisa, M.
(2007). KAAS: an automatic genome annotation and pathway reconstruction
server. \emph{Nucleic Acids Research} , \emph{35} (Web Server issue),
W182--W185. https://doi.org/10.1093/nar/gkm321

Mun, J.-H., Song, Y.-H., Heong, K.-L. \& Roderick, G.-K. (1999). Genetic
variation among Asian populations of rice planthoppers,\emph{Nilaparvata
lugens} and \emph{Sogatella furcifera} (Hemiptera: Delphacidae):
mitochondrial DNA sequences. \emph{Bulletin of Entomological Research} ,
\emph{89} , 245--253. https://doi.org/10.1017/S000748539900036X

Nam, H. Y., Kim, K. S., \& Lee, J. H. (2019). Population genetic
structure and putative migration pathway of \emph{Sogatella
furcifera}(Horvath)(Hemiptera, Delphacidae) in Asia. \emph{Bulletin of
entomological research} , \emph{109} (4), 453--462.
https://doi.org/10.1017/S0007485318000755

Nei M., \& Kumar S. Molecular evolution and phylogenetics. Oxford:
Oxford University Press; 2000.

Nazzi, F., Brown, S. P., Annoscia, D., Del Piccolo, F., Di Prisco, G.,
Varricchio, P., \ldots{} \&Pennacchio, F. (2012). Synergistic
parasite-pathogen interactions mediated by host immunity can drive the
collapse of honeybee colonies. \emph{PLoS Pathogen} , \emph{8} (6),
e1002735. https://doi.org/10.1371/journal.ppat.1002735

Nielsen, R., Williamson, S., Kim, Y., Hubisz, M. J., Clark, A. G., \&
Bustamante, C. (2005). Genomic scans for selective sweeps using SNP
data. \emph{Genome research} , \emph{15} (11), 1566--1575.
http://doi.org/10.1101/gr.4252305

Pavlidis, P., \selectlanguage{polish}Ž\selectlanguage{english}ivkovi\selectlanguage{polish}ć, \selectlanguage{english}D., Stamatakis, A., \& Alachiotis, N. (2013).
SweeD: likelihood-based detection of selective sweeps in thousands of
genomes. \emph{Molecular Biology and Evolution} , \emph{30} , 2224--34.
https://doi.org/10.1093/molbev/mst112

Pickrell, J. K., \& Pritchard, J. K. (2012). Inference of Population
Splits and Mixtures from Genome-Wide Allele Frequency Data. \emph{PLoS
Genetics} , \emph{8} (11), e1002967.
https://doi.org/10.1371/journal.pgen.1002967

Ruan, J., \& Li, H. (2020). Fast and accurate long-read assembly with
wtdbg2. \emph{Nature Methods} , \emph{17} (2), 155--158.
https://doi.org/10.1038/s41592-019-0669-3

Senol Cali, D., Kim, J. S., Ghose, S., Alkan, C., \& Mutlu, O. (2019).
Nanopore sequencing technology and tools for genome assembly:
computational analysis of the current state, bottlenecks and future
directions. \emph{Briefings in bioinformatics} , \emph{20} (4),
1542--1559.

Shang, F., Ding, B. -Y., Zhang, Y. -T., Wu, J. -J., Pan, S. -T., \&
Wang, J. -J. (2021). Genome-wide analysis of long noncoding RNAs and
their association with wing development in \emph{Aphis
citricidus}(Hemiptera: Aphididae).~\emph{Insect Biochemistry and
Molecular Biology} ,~\emph{139} , 103666.
https://doi.org/10.1016/j.ibmb.2021.103666

Shen, J. -H. (2003). Management of the whitebacked
planthoppers,\emph{Sogatella furcifera} in China: a mini review.
\emph{Chinese Journal Rice Science} , \emph{17} , 7--22.

Stamatakis, A. (2014). RAxML version 8: a tool for phylogenetic analysis
and postanalysis of large phylogenies. \emph{Bioinformatics} , \emph{30}
, 1312--1313. https://doi.org/10.1093/bioinformatics/btu033

Stanke, M., \& Morgenstern, B. (2005). AUGUSTUS: a web server for gene
prediction in eukaryotes that allows user-defined
constraints.\emph{Nucleic Acids Research} , \emph{33} (Web Server
issue), W465--W467. https://doi.org/10.1093/nar/gki458

Sun, J.-T., Jiang, X.-Y., Wang, M.-M., \& Hong, X.-Y. (2014).
Development of microsatellite markers for, and a preliminary population
genetic analysis of the white-backed planthopper. \emph{Bulletin of
Entomological Research} , \emph{104} (6), 765--773.
https://doi.org/10.1017/S0007485314000613

Turner, R., Song, Y. -H., \& Uhm, K. B. (1999). Numerical model
simulations of brown planthopper Nilaparvata lugens and white-backed
planthopper \emph{Sogatella furcifera} (Hemiptera : Delphacidae)
migration. \emph{Bulletin of Entomological Research} , \emph{89} ,
557--568. https://doi.org/10.1017/S0007485399000711

Wang, B., Kumar, V., Olson, A., \& Ware, D. (2019). Reviving the
transcriptome studies: an insight into the emergence of single-molecule
transcriptome sequencing. \emph{Frontiers in genetics, 10} , 384.
doi.org/10.3389/fgene.2019.00384

Wang, L., Tang, N., Gao, X.-L., Chang, Z.-X., Zhang, L.-Q., Zhou, G.-H.,
\ldots{} \& Wu, Q.-F. (2017). Genome sequence of a rice pest, the
white-backed planthopper (\emph{Sogatella furcifera}
).\emph{Gigascience} , \emph{6} (1).
https://doi.org/10.1093/gigascience/giw004

Wang, L., Tang, N., Gao, X. -L., Guo, D.-Y., Chang, Z.-X., Fu, Y.-T.,
\ldots{} \& Wu, Q.-F. (2016). Understanding the immune system
architecture and transcriptome responses to southern rice black-streaked
dwarf virus in Sogatella furcifera. \emph{Scientific Reports} , \emph{6}
, 36254. https://doi.org/10.1038/srep36254

Xi, Z., Ramirez, J. L., \& Dimopoulos, G. (2008). The Aedes aegypti toll
pathway controls dengue virus infection. \emph{PLoS Pathogen} ,\emph{4}
(7), e1000098. https://doi.org/10.1371/journal.ppat.1000098

Xie, J.-N., Guo, J.-J., Jin, D.-C., \& Wang, X.-J. (2014). Genetic
diversity of \emph{Sogatella furcifera} (Hemiptera: Delphacidae) in
China detected by intersimple sequence repeats. \emph{Journal of Insect
Science} , \emph{14} (1), 233. https://doi.org/10.1093/jisesa/ieu095

Xie, T., Zheng, J.-F., Liu, S., Peng, C., Zhou, Y.-M., Yang, Q.-Y., \&
Zhang, H.-Y. (2015). De novo plant genome assembly based on chromatin
interactions: a case study of \emph{Arabidopsis thaliana}
.\emph{Molecular Plant} , \emph{8} (3), 489--492.
https://doi.org/10.1016/j.molp.2014.12.015

Xu J, Cherry S. Viruses and antiviral immunity in Drosophila.
2014.\emph{Developmental and Comparative Immunology} , \emph{42} ,
67--84. https://doi.org/10.1016/j.dci.2013.05.002

Xue, J., Zhou, X., Zhang, C.-X., Yu, L.-L., Fan, H.-W., Wang, Z.,
\ldots{} \& Cheng, J.-A. (2014). Genomes of the rice pest brown
planthopper and its endosymbionts reveal complex complementary
contributions for host adaptation. \emph{Genome Biology} , \emph{15}
(12), 521. https://doi.org/10.1186/s13059-014-0521-0

Yan, L., Wang, X., Liu, H., Tian, Y., Lian, J.-M., Yang, R.-J., \ldots{}
\& Sheng, J. (2015). The genome of dendrobium officinale illuminates the
biology of the important traditional chinese orchid herb.\emph{Molecular
Plant} , \emph{8} (6), 922--934.
https://doi.org/10.1016/j.molp.2014.12.011

Yang, N., Dong, Z.-K., Chen, A.-D., Yin, Y.-Q., Li, X.-Y., \& Chu, D.
(2020). Migration of \emph{Sogatella furcifera} between the Greater
Mekong Subregion and northern China revealed by mtDNA and SNP. \emph{BMC
evolutionary biology} , \emph{20} (1), 1--10.
https://doi.10.1186/s12862-020-01722-4

Yang, Z. -H. (2007). PAML 4: phylogenetic analysis by maximum
likelihood. \emph{Molecular Biology and Evolution} , \emph{24} ,
1586--1591. https://doi.org/10.1093/molbev/msm088

Zambon, R. A., Nandakumar, M., Vakharia, V. N., \& Wu, L.-P. (2005). The
Toll pathway is important for an antiviral response in\emph{Drosophila}
. \emph{Proceedings of the National Academy of Sciences of the United
States of America} , \emph{102} (20), 7257--7262.
https://doi.org/10.1073/pnas.0409181102

Zhang, C.-X., Brisson, J., Xu, H.-J., (2019). Molecular mechanisms of
wing polymorphism in insects. \emph{Annual Review of Entomology}
,\emph{64} , 297--314.
https://doi.org/10.1146/annurev-ento-011118-112448

Zhang, L.-Q., Xu, W., Gao, X.-L., Li, W.-J., Qi, S., Guo, D.-Y.,
\ldots{} \& Wu, Q.-F. (2020). lncRNA Sensing of a Viral Suppressor of
RNAi Activates Noncanonical Innate Immune Signaling in \emph{Drosophila}
. \emph{Cell Host Microbe} , \emph{27} (1), 115--128 e118.
https://doi.org/10.1016/j.chom.2019.12.006

Zhou G.-H., Zhang. S.-G., Zou S.-F., Xu Z.-W., \& Zhou Z.-Q. (2010).
Occurrence and damage analysis of a new rice dwarf disease caused by
Southern rice black-streaked dwarf virus. \emph{Plant protection}
,\emph{36} , 144--146.

Zhou, G.-H., Xu, D.-G., Xu, D.-G., \& Zhang, M.-X. (2013). Southern rice
black-streaked dwarf virus: a white-backed planthopper-transmitted
fijivirus threatening rice production in Asia. \emph{Frontiers in
Microbiology} , \emph{4} , 270.

Zhu, J.-J., Jiang, F., Wang, X.-H., Yang, P.-C., Bao, Y.-Y., Zhao, W.,
\ldots{} \& Cui, F. (2017). Genome sequence of the small brown
planthopper,\emph{Laodelphax striatellus} . \emph{Gigascience} ,
\emph{6} (12), gix109. https://doi.org/10.1093/gigascience/gix109

\textbf{Figure Captions:}

\textbf{Figure 1} Schematic representation of the genome assembly and
genomic characteristics of the white-backed planthopper, \emph{Sogatella
furcifera} . (\textbf{a} ) Sequencing scheme and genome assembly
strategy adopted for \emph{S. furcifera} genome sequencing. A strategy
combining long-read PacBio sequencing, short-read Illumina sequencing,
and Hi-C sequencing schemes was employed to generate a high-quality and
complete genome assembly for \emph{S. furcifera} , `Surf\_2022'.
(\textbf{b} ) Assessments of BUSCO completeness for genome assemblies
for \emph{S. furcifera} . (\textbf{c} ) Hi-C-based heatmaps of DNA
interactions of each chromosome in `Surf\_2022'. This heatmap was
generated by Juicebox software using in situ Hi-C data from female WBPH
(\textbf{d} ) Genome landscape of the white-backed planthopper genome.
From the outer to the inner panel, these genome tracks represent the
following: I. Chromosome ideograms; II. Density of protein-coding genes;
III. Density of G+C content; IV. Density of long interspersed nuclear
elements (LINEs); V. Density of long terminal repeat retrotransposons
(LTR-RTs); VI. Density of DNA transposons; and VII. Density of
genome-wide single nucleotide polymorphism (SNP) sites.

\textbf{Figure 2} Gene family evolutionary analysis of \emph{Sogatella
furcifera} . (\textbf{a} ) Phylogenetic tree and orthologs in the
white-backed planthopper and 16 other arthropod species. The
phylogenetic tree was inferred from 860 strict single-copy orthologous
genes among 16 arthropod species using RAxML software. All the nodes
have at least 90\% bootstrap support. Numbers on internal nodes indicate
divergence time. MRCA is the abbreviation of the most recent common
ancestor. The histograms were subdivided to represent different
categories of orthology: 1:1:1 indicates single-copy universal genes
across these species. N:N:N stands for universal multiple-copy genes
across these species. The genes that were specific to Hymenoptera,
Lepidoptera, Diptera, and Coleoptera are also shown in the figure in
corresponding colors. Lineage-specific genes were also displayed for
each species. Other denotes genes that do not belong to these groups in
these species. Expansion and contraction analysis of gene families was
performed using CAF\selectlanguage{ngerman}É 3. The numbers of expanded (red) and contracted
(blue) gene families are marked in the corresponding branch. (\textbf{b}
) KEGG pathway enrichment analysis results of rapidly expanded gene
families of \emph{Sogatella furcifera} . Gene families with\emph{P}
\textless{} 0.05 were defined as rapidly evolving families. Pathways
with \emph{P} \textless{} 0.05 were considered significantly enriched.
(\textbf{c} ) Gene families rapidly were co-expanded and co-contracted
in the three rice planthoppers, white-backed planthopper (WBPH), brown
planthopper (BPH), and small brown planthopper (SBPH).

\textbf{Figure 3} NF-\selectlanguage{greek}κ\selectlanguage{english}B-related antiviral immune pathways
in\emph{Sogatella furcifera} . (\textbf{a} ) Toll signaling pathway.
Sp\selectlanguage{ngerman}ätzle (an endogenous ligand) binds to Toll receptors and transduces
signals to Cactus-Dif (dorsal-related immunity factor) through the
Toll-induced signaling complex, which contains MyD88, Tube and Pelle.
Cactus is phosphorylated and cleaved from Dif, allowing Dif to
translocate into the nucleus and regulate the transcription of effector
genes. (\textbf{b} ) Imd signaling pathway. The receptor molecules
multimerize upon activation and transmit the signal to the adaptor
molecule Imd. Imd activation recruits dFADD (Fas-associated death-domain
protein), which normally recruits a caspase-DREDD; however, dFADD was
not detected in white-backed planthopper and most of the Hemiptera
species involved in this study. Activation of DREDD leads to
polyubiquitination of Imd. TAK1 binds to the polyubiquitin chain and is
responsible for the assembly and activation of the IKK complex (IKK\selectlanguage{greek}β \selectlanguage{english}and
IKK\selectlanguage{greek}γ\selectlanguage{english}). The IKK complex mediates the phosphorylation of Relish.
Phosphorylated Relish is subsequently cleaved by DREDD, leaving the
C-terminal ankyrin repeats in the cytosol. The N-terminal DNA binding
domain of Relish translocates to the nucleus and regulates antimicrobial
peptides such as attacin C and diptericin B. The point at which viruses
interact with each of these pathways is currently unresolved. These two
simplified diagrams were referred to and modified from (Kingsolver et
al., 2013) and (Xu and Cherry, 2014). (\textbf{c} ) and (\textbf{d} )
Histograms of SRBSDV-P9-2 expression representing the virus propagation
levels through qPCR tests after the knockdown of representative immune
genes.

\textbf{Figure 4} Population genomic analysis of \emph{Sogatella
furcifera} in 12 geographical regions. (\textbf{a} ) Phylogenetic tree
of 12 geographic populations constructed from genome-wide SNP data based
on Nei's genetic distances. Abbreviations for sampling sites are as
follows: VNTB: Thai Binh, VNBG: Bac Giang, VNND: Nam Dinh, GXYZ:
Guangxi, GDGZ: Guangdongguangzhou, GDHF: Guangdonghaifeng, YNSD:
Yunnanshidian, YNMS: Yunnanmangshi, YNCJ: Yunnanchengjiang, HNBT:
Hainanbaoting, HNLS: Hainanlingshui, HNQH: Hainanqionghai. (\textbf{b} )
Minimum spanning network representing relationships among multilocus
genotypes (MLGs) of \emph{S. furcifera} identified from 12
geographically distinct populations. Each node represents a different
MLG. The size of a node represents the number of isolates in the same
MLG. Node colors represent population membership in populations.
(\textbf{c} ) Gene flow analysis across different geographic populations
in China and Vietnam. The possible gene flow events across geographic
populations were inferred using TreeMix when assuming the occurrence of
nine gene flow events. (\textbf{d} ) Discriminant analysis of principal
components (DAPC) analysis of all \emph{S. furcifera} isolates across 12
geographic populations. (\textbf{e} ) Genotype composition plot of
each\emph{S. furcifera} isolate in the 12 geographic populations from
DAPC analysis. The isolates are organized by geographic regions given at
the bottom.

\textbf{Figure 5} Genome-wide scan signals of selective sweeps in
white-backed planthopper (WBPH) populations. (\textbf{a} ) A plot of the
composite likelihood ratio (CLR) across the \emph{S. furcifera} genome.
The solid red line indicates the cutoff value for the 1\% outlier
regions with significant deviations from neutrality, indicating evidence
of recent selective sweeps. The genes with the strongest selective sweep
signals across the genome are labeled with different colors: purple
indicates the genes involved in morphogenesis, metamorphosis and other
crucial processes during insect development, blue represents the genes
involved in wing development, and red indicates the genes associated
with immunity. (\textbf{b} ) KEGG enrichment (\emph{P} value \textless{}
0.05) analysis of selected genes. (\textbf{c} ) GO enrichment (\emph{P}
value \textless{} 0.05) analysis of selected genes.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-1/Figure-1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-2/Figure-2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-3/Figure-3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-4/Figure-4}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-5/Figure-5}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

