Bulked Segregant RNA Sequencing Revealed Difference Between Virulent and Avirulent Brown Planthoppers

The brown planthopper (Nilaparvata lugens Stål, BPH) is one of the most devastating insect pests of rice (Oryza sativa L.), but BPH populations have varying degrees of virulence to rice varieties carrying different resistance genes. To help efforts to characterize these variations we applied bulked segregant RNA sequencing (BSR-seq) to identify differentially expressed genes (DEGs) and genetic loci associated with BPH virulence to YHY15 rice plants carrying the resistance gene Bph15. BPHs that are highly virulent or avirulent to these plants were selected from an F2 population to form two contrasting bulks, and BSR-seq identified 751 DEGs between the bulks. Genes associated with carbohydrate, amino acid and nucleotide metabolism, the endocrine system, and signal transduction were upregulated in the avirulent insects when they fed on these plants. The results also indicated that shifts in lipid metabolism and digestive system pathways were crucial for the virulent BPHs’ adaptation to the resistant rice. We identified 24 single-nucleotide polymorphisms (SNPs) in 21 genes linked with BPH virulence. Possible roles of genes apparently linked to BPH virulence are discussed. Our results provide potentially valuable information for further studies of BPH virulence mechanisms and development of robust control strategies.

Introduction

The brown planthopper (Nilaparvata lugens Stål, BPH) is a monophagous piercing-sucking herbivore that has become the most destructive insect pest of rice (Oryza sativa L.) (Sogawa, 1982; Hibino, 1996; Guo et al., 2018). BPHs occur widely throughout Asia and can migrate several 100 km across national borders (Jing et al., 2017). When feeding on rice plants, BPH nymphs and adults ingest phloem sap by inserting a stylet into their vascular tissues, thereby directly damaging the plants (Sogawa, 1982). They can also directly damage the plants by transmitting pathogenic viruses (Jena and Kim, 2010; Guo et al., 2018). Rice varieties that are resistant to BPH can inhibit the pest’s feeding, growth, and fecundity (Zheng et al., 2021). In recent decades, many BPH resistance genes and QTLs have been detected in cultivated rice and wild species (Guo et al., 2018; Shi et al., 2021; Zhou et al., 2021). However, new BPH biotypes may evolve that are virulent to resistant rice varieties (Claridge and Den Hollander, 1980; Saxena et al., 1991). A virulent biotype is defined as a BPH population that can damage a rice variety that was previously resistant (Claridge and Den Hollander, 1980; Jena and Kim, 2010). The first BPH-resistant variety (IR26 carrying the Bph1 resistance gene) was developed by the International Rice Research Institute (IRRI) in 1973. However, within several years BPH biotype 2 had adapted and become virulent to IR26 in the field (Jena and Kim, 2010). Rice varieties resistant to biotype 2, carrying Bph2, were subsequently widely deployed in the field, but BPH biotype 3, which can overcome Bph2-mediated resistance, was detected in 1981 (Pathak and Heinrichs, 1982; Jena and Kim, 2010). Several additional BPH biotypes have been developed by rearing biotype 1 insects on specific resistant rice varieties in the laboratory (Claridge and Den Hollander, 1982). Analyses of various biotypes have clearly shown that BPH’s virulence to resistant rice varieties is a genetic trait (Kobayashi, 2016), but the mechanisms underlying BPH virulence are poorly understood.

Previous studies have shown that BPHs exhibit distinct responses to feeding on resistant plants carrying different resistance genes (Kobayashi, 2016; Yue et al., 2019). For instance, their responses to the resistant rice variety Mudgo, which carries the Bph1 resistance gene, include activation of various metabolic, digestive, and absorption pathways (Ji et al., 2013). Wang et al. (2015) also found that feeding on the resistant rice variety B5 induced significant changes in expression of genes involved in metabolism of lipids, amino acids, carbohydrates, nucleotides, and TCA cycle intermediates in the salivary glands of BPH insects. Zhang et al. (2019) showed that expression levels of genes related to detoxification and starvation responses were significantly upregulated in BPHs that fed on resistant Bph6-transgenic rice. In addition, Zheng et al. (2020) observed stronger expression of genes involved in lipid metabolism in BPHs fed on NIL-Bph6 and NIL-Bph9 resistant rice plants than in insects fed on susceptible wild type rice. These studies have provided valuable insights. However, most data on BPH responses have been obtained by temporarily feeding avirulent insects on different rice varieties, and there have been few comparisons of virulent and avirulent BPHs fed on resistant varieties.

Various techniques have been used to study plant-insect interactions (Zogli et al., 2020). For instance, genome sequencing of the pea aphid showed that most chemosensory genes had undergone rapid expansion in the genome under positive selection (Smadja et al., 2009). Genome sequencing has detected the expansions of detoxification gene families in Tetranychus urticae (Grbić et al., 2011). The genome-wide analysis found the expansion of Cytochrome P450 (CYP) genes in the rice striped stem borer in response to rice defense secondary metabolites (Wang et al., 2014). MS/MS de novo peptide sequencing has found several enzymes responsible for digesting carbohydrates, proteins, and lipids in the midgut lumen of the cotton bollworm (Pauchet et al., 2008). In addition, insects have evolved effectors to adapt to host plants (Chen et al., 2019). The functional genomics approach has found effectors required for thrips to reproduce on tomatoes (Abd-El-Haliem et al., 2018). Additionally, some insect avirulence genes have been identified in Hessian fly associated with virulence to wheat R genes H6 and H13 by bulked segregant analysis (BSA)-seq analysis (Navarro-Escalante et al., 2020). A useful approach for studying disease resistance in plants is BSA, initially developed by Michelmore et al. (1991). BSA is based on the theory that different genotypes are likely to have differing phenotypic extremes, which can be used to rapidly identify genetic markers linked to genomic regions related to targeted traits (Michelmore et al., 1991; Wang et al., 2013). However, due to complications caused by repetitive genomic sequences and high genome sequencing depth requirements, BSA is not suitable for analyses of populations with numerous mutants and large genomes (Liu et al., 2012). Another approach, RNA sequencing (RNA-seq), is widely applied in next-generation sequencing (NGS) efforts to quantify relative quantities of specific transcripts from read counts (Górczak et al., 2020). RNA-seq reads can also be used to identify DNA sequence polymorphisms, such as single nucleotide polymorphisms (Chepelev et al., 2009). Bulked segregant RNA-seq (BSR-seq) couples BSA with RNA-seq (Liu et al., 2012) and can be used not only to identify differentially expressed genes (DEGs) but also single-nucleotide polymorphisms (SNPs) that are tightly linked to targeted traits (Du et al., 2017). BSR-seq has been successfully applied in analyses of many species. For instance, BSR-seq has identified 1,255 DEGs between catfish that are resistant and susceptible to enteric septicemia, and 56,419 SNPs of 4,304 unique genes of the fish (Wang et al., 2013). In addition, it has identified 18 high-probability SNPs and six genes that may be involved in responses to waterlogging stress in maize (Du et al., 2017), as well as the maize glossy13 gene, which encodes a putative ABC transporter and is crucial for accumulation of epicuticular waxes (Li et al., 2013). As it combines advantages of BSA and RNA-seq it is highly efficient for analyses of complex traits involving multiple genes (Du et al., 2017; Li et al., 2020).

The BPH resistance gene Bph15 was introduced into O. sativa from the wild rice species Oryza officinalis (Yang et al., 2004; Lv et al., 2014). Rice variety YHY15 is a recombinant inbred line (RIL) selected from the RI93 × TN1 F2 population and carries a single resistance gene Bph15 (Yang et al., 2004). BPH population, designed biotype Y, was developed by rearing biotype 1 BPHs on YHY15 from January 2007 (Jing et al., 2012). Biotype Y can overcome the resistance conferred by Bph15, grow and reproduce normally on YHY15 plants. Thus, this combination of insects and plants provides a suitable model system for examining virulence mechanisms, and in this study we applied BSR-seq to identify DEGs and SNPs associated with BPH virulence to rice. Comparison of virulent and avirulent BPH extreme bulks after feeding on resistant YHY15 rice plants for 48 h identified 751 DEGs, with illuminating differences in transcriptomic profiles. Moreover, we identified 24 polymorphic SNPs related to BPH virulence located on six chromosomes, and 21 genes that may play key roles in virulence mechanisms. The results may facilitate further study of the mechanisms that enable BPH to overcome host resistance.

Materials and Methods

Insect and Plant Materials

The brown planthopper biotype 1 insects were reared on the susceptible rice variety Taichung Native 1 (TN1), which does not contain any BPH resistance gene. BPH biotype Y was developed by forcing biotype 1 BPHs to feed and reproduce on the resistant rice variety YHY15 for multiple generations (Jing et al., 2012). YHY15 is a RIL carrying the Bph15 resistance gene (Yang et al., 2004). The BPHs used in the experiments reported here were reared on 1-month-old rice seedlings growing in 24 × 19 cm plastic buckets in controlled environmental conditions (26 ± 1°C, 16 h light/8 h dark cycles) at Wuhan University Institute of Genetics.

Brown Planthopper Weight Gain and Honeydew Excretion Assay

The brown planthopper weight gain and honeydew excretion were measured as previously described (Shi et al., 2021). Newly emerged female adults were weighed using an AUW120D electronic balance (AUW120D, Shimadzu, Japan) then introduced into separate pre-weighed parafilm sachet (2 × 2.5 cm) fixed to leaf sheaths of a 4-week-old rice plant. After 48 h the insects were carefully removed from the sachets, and then the insect and the honeydew in each sachet were weighed. Each insect’s weight gain was determined from its weight before and after ingestion, and its weight gain ratio was calculated by dividing its weight gain by its initial weight. Both BPH weight gain and honeydew excretion assays were repeated three times with 25 BPHs for each replicate.

Brown Planthopper Survival Rate

The brown planthopper survival rates were determined as previously described (Shi et al., 2021). Briefly, rice plants were grown in plastic pots (12 cm tall, 9 cm diameter, one plant per pot). When they were 1-month-old the pots were covered with plastic cages and 10 third-instar nymphs were released into each cage to infest the plants. The number of surviving BPH nymphs in each cage was recorded every day until 10 days after their introduction, and survival rates were calculated by dividing the number of surviving nymphs in each cage by the number initially released. The experiment was repeated three times, with 20 BPHs per replicate.

Formation of Brown Planthopper Extreme Bulks

An extremely virulent BPH female biotype Y adult, which gained 2.06 mg weight in the weight gain assay described above, was mated with an avirulent male adult of biotype 1. The F1 offspring were allowed to intermate and female F2 adults resulting from their mating, with similar body sizes, were subjected to the assay. Forty that gained 1.89–2.8 mg weight were classified as virulent to YHY15, and 40 that lost 0.4–1.04 mg were classified as avirulent. These sets were used to form a Bph15-virulent bulk (vB15) and a Bph15-avirulent bulk (aB15), respectively, which were used to prepare bulked RNA samples for the BSR-seq analysis.

Illumina cDNA Library Construction and RNA Sequencing

Each BPH individual in the extreme bulks was dissected into two parts under an Olympus stereoscopic microscope, the heads with intestines and the remaining carcasses. The two resulting bulks from the sets of 40 virulent and avirulent F2 individuals were each randomly divided into four replicates, with parts of 10 insects in each replicate. The tissue samples in each replicate were mixed, and total RNA was extracted from them using a mirVana miRNA Isolation Kit (Ambion-1561, Austin, TX, United States). After testing its integrity using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States), sequencing libraries were constructed using a TruSeq Stranded mRNA LTSample Prep Kit (Illumina, San Diego, CA, United States) following the manufacturer’s instructions. The cDNA library preparations were then sequenced using an Illumina HiSeq 4000 sequencing platform. Raw reads were processed using Trimmomatic (Pérez-Rubio et al., 2019). To obtain clean data, reads containing poly-n and low-quality reads were removed. Then the clean reads were used for the BSR-seq analysis.

SNP Calling and Bulked Segregant RNA-Seq Analysis

The reads were mapped to the latest chromosome-level BPH reference genome (BioProject accession no. PRJNA591478) (Ma et al., 2021) using Burrows–Wheeler Aligner (BWA, v. 0.7.5a) (Rajan-Babu et al., 2021). The bwa index command and BWA–MEM were used to create reference genome index and map the reference genome, respectively. Then, the mapped reads were sorted and indexed by SAMtools (v. 1.3.1) (Danecek et al., 2021), following recommended procedures.¹ SNPs were called using SAMtools and BCFtools (v. 1.3.1) (Lam et al., 2020). The Variant Call Format (VCF) files are shown in Supplementary Materials. Then SnpEff v. 4.1g (Cingolani et al., 2012) was used to annotate and predict the effects of genetic variants on genes. And SNP index was calculated at each SNP position to indicate the proportion of reads harboring SNPs that differed from the reference sequences (Abe et al., 2012). SNP positions with SNP index less than 0.3 were eliminated, and those with depth less than four were excluded, as these SNPs may be due to sequencing or alignment errors (Hisano et al., 2017; Feng et al., 2020; Li et al., 2020). After filtration, ΔSNP-index value (Takagi et al., 2013; Li et al., 2020) was calculated for each position by subtracting the SNP index of the aB15-bulk from that of the vB15-bulk. The average ΔSNP-index was calculated using a sliding window analysis with a window size of 3-Mb and a slide size of 1-Mb. We computed confidence intervals of ΔSNP-index values for all SNP positions under the null hypothesis that no QTLs were present.

Differential Expression Gene Analysis

The clean reads were mapped to the BPH reference genome (Ma et al., 2021) using Hisat2 (v. 2.2.1.0) software² (Kim et al., 2015). To quantify gene expression levels, fragments per kilobase of transcript per million mapped reads (FPKM) values (Roberts et al., 2011) were calculated using Cufflinks (Trapnell et al., 2010), and read counts of each gene were determined by htseq-count (Srinivasan et al., 2020). DEGs were clustered using the DESeq package (Anders and Huber, 2012) with thresholds of log₂FC ≥ 1 and P < 0.05 for significantly differential expression. The DEGs’ expression patterns were explored by hierarchical cluster analysis and they were functionally annotated by comparing with entries in the Nr, eukaryotic orthologous groups (KOG), Swiss-Prot, Gene ontology (GO), Pfam, and Kyoto Encyclopedia of Genes and Genomics (KEGG) databases (with P ≤ 0.05 as the threshold in GO term and KEGG pathway enrichment analyses).

Real-Time Quantitative PCR Analysis

For expression analyses of candidate genes in Bph15-virulent (Vr-B15) and Bph15-avirulent (Avr-B15) BPHs. The weight gain of newly emerged female adults from biotype Y was measured after feeding on YHY15 for 48 h. The virulent insects (average weight gain in 2.0 mg) and the avirulent (weight gain –0.4 mg) were used for qRT-PCR, respectively. Total RNA was isolated from the insects using a RNAiso Plus kit (Takara, Kyoto, Japan) following the manufacturer’s instructions. RNA samples were reverse-transcribed into cDNA using a PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Kyoto, Japan). The resulting sequences were then used in qRT-PCR analysis with a CFX96 Touch™ Real-Time System and iTaq Universal SYBR Green Supermix Kit (Bio-Rad, Hercules, CA, United States) according to the manufacturer’s protocols. Relative expression levels were calculated using the 2^–ΔΔCt method with the BPH housekeeping gene actin serving as an internal control. Three independent biological replicates were employed in each experiment.

Statistical Analysis

The significance of differences in recorded variables was assessed using Student’s t-test or Tukey’s honestly Significant Difference test implemented in IBM SPSS Statistics v. 23.0. Differences were considered significant, very significant and highly significant (indicated by *, ^**, ^*** in the figures) if P < 0.05, P < 0.01, and P < 0.001, respectively.

Results

The Performance of Brown Planthopper Insects on YHY15 Rice Plants

To verify the BPH resistance conferred by Bph15, YHY15 plants carrying Bph15 were used for seedling stage insect resistance analysis, with susceptible TN1 plants as controls. After infestation by biotype 1 BPHs for 7 days, the TN1 plants were dead, but the YHY15 plants survived well, proving YHY15’s resistance to biotype 1 BPHs (Figures 1A,B and Supplementary Figure 1). Moreover, the mortality rate of biotype 1 insects was also significantly higher on YHY15 (24.44%) than on TN1 (7.22%) plants (Figure 1C). These results clearly indicate that the YHY15 plants were resistant to BPH.

FIGURE 1

Figure 1. Bph15-mediated resistance verification and virulence assessment of BPH biotype Y. (A) Photo showing the difference in TN1 and YHY15 rice plants’ BPH resistance. (B) BPH-resistance scores of TN1 and YHY15 plants (negatively related to BPH resistance) obtained from observations 7 days after BPH infestation (n = 3 independent experiments, each with 15 rice plants). (C) Mortality of BPH insects fed on TN1 and YHY15 plants after 48 h. Weight gain (D) and growth rate (E) of BPH biotypes 1 and biotype Y feeding on YHY15 rice plants for 48 h. Three independent experiments were repeated with 25 BPHs per replicate. (F) Amounts of honeydew excreted by newly emerged female adults of biotype 1 and biotype Y on YHY15 plants after 48 h. Three independent experiments were performed with 25 BPHs per replicate. All data presented are means ± SE of three replicates. In the box charts, the bottom and top edges of the box charts represent the 25th and 75th percentiles, respectively, and center values are medians. The whiskers mark ranges of the data, excluding outliers. **, *** on the bars indicate significance at P < 0.01, and P < 0.001 (according to Student’s t-test), respectively.

To further investigate the virulence of BPH biotype Y to YHY15 rice, we observed the performance of adult female biotype 1 and Y insects (with similar initial body weight, ranging from 1.8 to 2.7 mg) on YHY15 plants. After feeding on YHY15 plants for 48 h, the biotype Y and 1 insects had average weight gain of 1.08 and 0.03 mg, respectively (Figure 1D), with average growth rate of 51.58 and 2.24%, respectively (Figure 1E). In addition, the biotype Y and 1 BPHs had average honeydew excretion of 42.78 and 3.55 mg, respectively (Figure 1F). These results indicate that the biotype Y has adapted to the resistant YHY15 rice plants.

Construction of Extremely Virulent and Avirulent Brown Planthopper Bulks

To construct the virulent and avirulent BPH bulks, 504 BPH female adults of biotype Y were randomly selected and we recorded the weight gain of each insect after feeding on YNY15 plants for 48 h. Their weight gain of BPH insects was continuously distributed in range of –0.79–2.06 mg (Figure 2A), indicating significant heterogeneity in virulence in the biotype Y population. Therefore, a virulent female with 2.06 mg weight gain from biotype Y (indicated by the red asterisk in Figure 2A) was selected and crossed with an avirulent male of the biotype 1 population. The emerging male and female offspring were then intercrossed to breed a F2 generation for bulk construction.

FIGURE 2

Figure 2. Construction of the two extreme BPH bulks. (A) Frequency distribution of weight gain of 504 biotype Y BPHs after feeding on YHY15 plants for 48 h. A single extremely virulent BPH female, with 2.06 mg weight gain (indicated by the red asterisk), was selected. (B) Frequency distribution of weight gain of F2 BPHs (n = 1202) fed on YHY15 plants for 48 h. The histogram shows the BPHs’ weight gain. The chart above it shows the range of BPH weight gain, with 25th and 75th percentiles (excluding outliers) indicated by right and left edges of the box, respectively, and the median within the box. Panels (C–E) show the weight gain, weight gain ratio, and honeydew excretion, of the two bulks after feeding on YHY15 for 48 h (three independent experiments were performed). vB15-bulk and aB15-bulk refer to BPHs extremely virulent and extremely avirulent to YHY15, respectively. ***P < 0.001 (Student’s t-test).

In total, the virulence of 1202 newly emerged F2 female adults was assessed by allowing them to feed on resistant YHY15 plants for 48 h and recording their weight gain, which ranged from –1.04 to 2.8 mg (Figure 2B), clearly showing that their virulence widely varied. We selected 40 extremely virulent and 40 avirulent BPHs with weight gains ranging from 1.89 to 2.8 mg and weight losses ranging from 0.4 to 1.04 mg, respectively (Figure 2C), and average weight gain ratios of 96.49 and –22.27%, respectively (Figure 2D). In addition, BPH honeydew excretion is a reliable indicator of BPH feeding activity. It has been applied to determine the virulence of the biotypes (Claridge and Den Hollander, 1980). The virulent BPHs produced 25.62–94.45 mg honeydew, and the avirulent 0–0.77 mg honeydew, respectively (Figure 2E). Hence, the 40 extremely virulent F2 insects and 40 extremely avirulent counterparts were designated the Bph15-virulent (vB15) bulk and Bph15-avirulent (aB15) bulk, respectively (Supplementary Figure 2).

Transcriptome Annotation of the Assembled Unigenes

We subjected the two extreme BPH bulks to deep RNA-sequencing. For this, the 40 F2 individuals in each extreme bulk were randomly divided into four replicates with 10 insects in each replicate, and each insect was divided into two parts (heads + intestines and remaining carcasses). RNA transcripts were extracted from the pooled pairs of tissue samples of each set of 10 BPHs then sequenced using an Illumina HiSeq 4000 sequencing platform, generating eight transcriptome libraries with 59.6 G raw data in total. After filtering out low-quality sequences, incorrect adapters, and redundant sequences, 55.67 G clean data were obtained. In addition, 82.09–86.03% of clean reads were mapped to the latest chromosomal-level BPH reference genome (Ma et al., 2021), and the Q30 percentage exceeded 91.58% (Table 1).

TABLE 1

Table 1. Summary of the RNA sequencing data.

Results of the correlation analysis showed that the gene expression profiles of the four biological replicates of each bulked RNA sample correlated well (Supplementary Figure 3 and Supplementary Table 1). Principle component analysis (PCA) showed that the first two principal components (PC1 and 2) explained 61 and 12.7% of the variance, respectively. Samples within each group clustered together (Supplementary Figure 4). In particular, 18,559 unigenes were identified based on the clean reads obtained from the two BPH extreme bulks. Based on the sequence homology, 10,450 unigenes were assigned to KOG functional classifications. These included 25 KOG categories, most commonly general function prediction only (21.48%), followed by: signal transduction mechanisms (14.08%); post-translational modification, protein turnover, chaperones (8.871%); transcription (6.134%); intracellular trafficking, secretion, and vesicular transport (5.33%); lipid transport and metabolism (4.813%); amino acid transport and metabolism (4.584%); and translation, ribosomal structure and biogenesis (3.962%) (Figure 3). Transcripts of numerous genes involved in basic processes such as translation and transcription, as well as functions including lipid transport, signal transduction, and post-translational modification mechanisms were detected.

FIGURE 3

Figure 3. Histogram of numbers of genes in orthologous group classes. In total 10,450 unigenes were annotated into 25 KOG categories and the number of unigenes in each specific functional clump is indicated on the y-axis.

Transcriptomic Profiles of Avirulent and Virulent Brown Planthoppers

Transcriptomic profiles of the vB15- and aB15-bulks were compared to identify DEGs with fold-change (FC) ≥ 2 and P < 0.05 (Supplementary Figure 5). In total, 751 DEGs were detected (Supplementary Table 2), of which 418 were more strongly expressed in the aB15-bulk and 333 in the vB15-bulk (Figure 4A). These findings indicate that feeding on the YHY15 plants induced stronger transcriptomic changes in the avirulent BPHs than in the virulent BPHs. We subsequently subjected to the DEGs to gene ontology (GO) analysis to identify the key processes that were enriched in each extreme bulk. The GO analysis (level 2) illustrated that GO terms that were mainly enriched in both aB15 and vB15-bulks were involved in single-organism, cellular, metabolic, binding, and catalytic activity processes (Supplementary Figure 6 and Supplementary Table 3). Additionally, KEGG analysis (level 2) revealed that the 418 DEGs more strongly expressed in the aB15-bulk were mainly involved in: carbohydrate, amino acid and nucleotide metabolism; signal transduction; the endocrine system; and protein folding, sorting and degradation (Figure 4B). These changes indicate that Bph15-mediated resistance mechanisms induced substantial changes in diverse physiological processes of avirulent BPHs. The 333 DEGs more strongly expressed in the vB15-bulk were mainly assigned to lipid metabolism and digestive system categories (Figure 4B), suggesting that these pathways are crucial for BPH adaption to Bph15-mediated resistance.

FIGURE 4

Figure 4. Functional analyses of the differentially expressed genes (DEGs). (A) Volcano plots of DEGs (P < 0.05 and log2FC = 1.) between the vB15- and aB15-bulks, based on log2FC and FDR values. Green and red scatterpoints indicate DEGs more strongly expressed in the aB15-bulk and vB15-bulk, respectively. The results are based on comparison of expression of genes in vB15-bulk and aB15-bulk. (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) classification of the DEGs more strongly expressed in the aB15- and vB15-bulks, respectively. The DEGs were assigned to KEGG pathways associated with: (A) Cellular processes; (B) Environmental information processing; (C) Genetic information processing; (D) Human diseases; (E) Metabolism; (F) Organismal systems. The DEGs more strongly expressed in the aB15-bulk were assigned to 42 KEGG categories, and those more strongly expressed in the vB15-bulk were assigned to 36 categories. (C,D) The 20 most highly enriched KEGG pathways in the aB15 and vB15 bulks, respectively.

We also subjected the DEGs to KEGG enrichment analysis (with a P < 0.05 threshold) to further pinpoint key pathways in each bulk. The 418 more strongly expressed in the aB15-bulk were mainly involved in: purine, carbohydrate, and amino acid metabolism; longevity regulation; and MAPK signaling pathways. Genes associated with six amino acids (Gly, Ser, Thr, Val, Leu, and Ile) were more strongly expressed in the aB15-bulk after feeding on resistant YHY15 rice plants. Genes of the glycolysis/gluconeogenesis, citrate cycle (TCA cycle), and pentose phosphate pathways of carbohydrate metabolism categories were significantly enriched in the avirulent BPHs (Figure 4C and Supplementary Table 4). In contrast, assignations of the 333 DEGs more strongly expressed in the vB15-bulk indicated significant enrichment of lipid metabolism (including steroid hormone and unsaturated fatty acid biosynthesis, fatty acid elongation, arachidonic acid metabolism, and glycerophospholipid metabolism pathways) in the virulent insects. Pancreatic secretion, vitamin digestion and absorption, and protein digestion and absorption pathways of the digestive system category were also significantly enriched in the vB15-bulk (Figure 4D and Supplementary Table 4). These results indicate that the virulent BPHs’ adaption to YHY15 is related to lipid metabolism, protein digestion, and absorption pathways.

Verification of Differentially Expressed Genes in Key Pathways

To verify the reliability of the RNA sequencing data, we validated 30 randomly selected DEGs by quantitative reverse transcription-PCR (qRT-PCR) analysis with gene-specific primers (Supplementary Table 5). Expression levels of these 30 randomly selected DEGs (13 and 17 that were more strongly expressed in the vB15-bulk and aB15-bulk, respectively) were determined in both bulks, normalized with respect to actin expression, and then compared with the corresponding transcript read counts. Their expression levels according to the RNA-seq and qRT-PCR analyses correlated extremely well, with an R² value of 0.9574 (Figure 5A), clearly confirming the credibility of our RNA-seq data.

FIGURE 5

Figure 5. Validation of the transcriptome sequencing quality. (A) Correlation of the RNA-seq and RT-qPCR results. Expression levels of 30 randomly selected DEGs, including 13 and 17 more strongly expressed in the vB15-bulk and aB15-bulks, according to analyses of the RNA samples (normalized with respect to the BPH β-actin gene; three biological replicates). Fold-change values of these genes obtained from the RNA-seq and RT-qPCR analyses are strongly positively correlated. (B) Expression levels of nine DEGs related to lipid metabolism in Avr-B15 and Vr-B15 (avirulent and virulent BPHs with 0.4 mg weight loss and 2.0 mg weight gain after feeding on a YHY15 plant for 48 h, respectively). Presented data are means ± SE of three replicates. Error bars indicate SEM, n = 3 independent experiments. Asterisks above the bars indicate significant differences, according to Student’s t-test (*P < 0.05; **P < 0.01; ***P < 0.001).

We also analyzed nine DEGs in the lipid metabolism pathway that were significantly enriched using qRT-PCR (Supplementary Table 5). These nine DEGs’ annotations indicate that they are UDP-glucuronosyltransferases (Supplementary Table 2), which play diverse roles, including detoxification (Bock, 2016), by catalyzing the transfer of sugar molecules to lipid and protein acceptors (Zielinski et al., 2016). Expression levels of these genes were examined, by qRT-PCR, in representatives of the Bph15-virulent (Vr-B15) and Bph15-avirulent (Avr-B15) insects of biotype Y, with 2.0 mg weight gain and 0.4 mg weight loss after feeding on YHY15 plants for 48 h, respectively. The results showed that all these genes were expressed significantly more strongly in the Vr-B15 insect than the Avr-B15 insect (Figure 5B). These results are consistent with the pathway enrichment results, indicating that lipid metabolism pathways were regulated differentially in the virulent and avirulent insects in response to rice plants carrying the Bph15 resistance gene.

Identification of Single-Nucleotide Polymorphisms Linked to Virulence

In total, 63127 SNPs with high confidence were obtained from the transcriptomic data of the two bulks. Application of the ΔSNP-index method identified 24 SNPs localized in 21 genes on six chromosomes at a significance level of 99% (Supplementary Table 6; Feng et al., 2020; Li et al., 2020; He et al., 2022). Moreover, two chromosome regions harboring SNPs were identified on chromosome 3 (73413752-83702207 and 88874910-94058321), as shown in Figure 6 and Supplementary Figure 7. In addition, 33 of the 751 DEGs were located near identified SNPs linked to virulence (Supplementary Table 7). These results clearly suggest that multiple loci control BPH virulence to plants with Bph15-mediated resistance.

FIGURE 6

Figure 6. Distribution of polymorphic SNPs on chromosomes. The horizontal lines indicate positions of SNPs with a significance level of 99%. The green blocks indicate chromosome regions linked to BPH virulence.

Putative functions of the 21 genes hosting the 24 identified SNPs were annotated (Supplementary Table 8). One of these genes, LOC111057305, encodes a fatty acyl-CoA reductase (FAR), which plays a crucial role in insect pheromone biosynthesis (Teerawanichpan et al., 2010). Another, LOC111052455, encodes a member of the aldo–keto reductase 6 (AKR6) family involved in channel trafficking and membrane localization (Barski et al., 2009). One of these genes, LOC111054236, encodes a dipeptidyl peptidase 4 (DPP 4) protein, a vital digestive peptidase in Tenebrio molitor larvae (Tereshchenkova et al., 2016). Another, LOC111046750, encodes a peroxidasin homolog protein involved in basement membrane biogenesis, tissue development, and innate immune defense (Nelson et al., 1994). And LOC120354300 encodes a zinc finger protein, which has diverse functions, including lipid binding, transcriptional activation, regulation of apoptosis, and protein folding and assembly (Laity et al., 2001). Others include LOC111055290, an integrator complex subunit implicated in adipose differentiation (Otani et al., 2013), and LOC111048257, encoding a basement membrane specific heparan sulfate proteoglycan core protein. These candidate genes encode enzymes, zinc finger proteins, and other functional proteins that may participate in BPH virulence to plants with Bph15-mediated resistance.

Discussion

Previous studies have revealed significant genetic variations in virulence within BPH populations (Den Hollander and Pathak, 1981; Tanaka, 1999). Accordingly, this study found substantial variation in the virulence of randomly selected individuals of BPH biotype Y (Figure 2A). Thus, direct comparison of wild BPH populations may not readily identify genetic factors of BPH virulence (Via, 1990; Kobayashi, 2016). The heterogeneity of the parental BPH population also reportedly complicates inheritance of BPH virulence (Sogawa, 1981). Therefore, we mated a virulent female and an avirulent male to produce a F1 generation, and F1s were intermated to produce the F2 generation. Two BPH extreme bulks were constructed from the F2 population by the BSR approach. Virulence assessment validated these two extreme bulks (Figure 2), then 751 DEGs and 24 SNPs in 21 candidate genes were identified by BSR-seq analysis, confirming its suitability for studying BPH biotype Y’s virulence to rice plants with Bph15-mediated resistance. To our knowledge, this is the first reported use of the BSR-seq method to identify candidate genes involved in BPH virulence to rice resistance.

Functional analyses of the identified DEGs showed that feeding on plants with Bph15-mediated resistance induced more extensive transcriptomic changes in the avirulent BPHs than in the virulent insects. Carbohydrates are the main chemical components in the phloem sap of rice and essential nutrients for phloem-sucking insects (Foyer et al., 2007). In this study, we detected 42 DEGs associated with carbohydrate metabolism in the two BPH bulks. More of these were expressed more strongly in avirulent BPHs than in virulent BPHs (Figure 4B). The expression patterns also showed that the glycolysis/gluconeogenesis, citrate cycle (TCA cycle), and pentose phosphate pathways were significantly enriched in the avirulent BPHs after feeding on resistant YHY15 rice plants (Figure 4B). These findings indicate that feeding on YHY15 plants induces upregulation of carbohydrate metabolism in avirulent BPH. We also found that it induced significant increases in expression of genes involved in synthesis of six amino acids (Gly, Ser, Thr, Val, Leu, and Ile) in the avirulent BPHs (Figure 4C). This has clear physiological implications as amino acids are the raw materials for protein synthesis and participate in myriads of fundamental processes in insects, including energy provision (Hargrove, 1976), sclerotization of newly formed cuticle (Andersen, 2007), and detoxification of harmful exogenous substances (Novoselov et al., 2015). Similarly, Peng et al. (2017) found that levels of six amino acids (Val, Leu, Ser, Thr, Pro, and Gln) were higher in BPHs fed on resistant YHY15 rice plants than in counterparts fed on susceptible TN1 rice. Moreover, Liu et al. (2017) found that levels of 19 amino acids were lower in BPH nymphs after feeding on YHY15 rice than after feeding on TN1 plants for 12 h. These results suggest that the higher amino acid levels may contribute to the feeding and survival of avirulent BPHs on resistant YHY15 rice. MAPK signaling, purine metabolism, and longevity regulating pathways were also enriched in the avirulent BPHs after feeding on the resistant rice (Figure 4C). Collectively, these findings suggest that feeding on rice with Bph15-mediated resistance may perturb diverse physiological processes in avirulent BPHs.

In contrast, feeding on the YHY15 plants more specifically altered lipid metabolism in virulent biotype Y BPHs, in accordance with previous findings that lipid mobilization in their bodies is crucial for BPHs feeding on resistant rice plants (Yue et al., 2019). In this study, we found that it induced stronger expression of lipid metabolism-related pathways in virulent BPHs than in avirulent BPHs. Biosynthesis of steroid hormones and unsaturated fatty acids, fatty acid elongation, arachidonic acid metabolism, and glycerophospholipid metabolism pathways were significantly enriched in the virulent-BPH bulk (Figure 4D). These results indicate that lipid metabolism pathways are important for BPH virulence to Bph15 resistance, and are consistent with previous reports that lipid metabolism in BPH fat bodies is activated when the insects feed on rice plants carrying Bph6 or Bph9 resistance genes (Zheng et al., 2020). In addition, efficient digestion and nutrient uptake are vital for insects to acquire essential nutrients from their diet (Holtof et al., 2019). Comparative analysis of salivary gland transcriptomes of two BPH populations derived from TN1 and Mudgo has revealed that genes related to “digestion and absorption” might be associated with BPH virulence (Ji et al., 2013). Results presented here indicate that feeding on the resistant plants induced stronger expression of pathways involved in the digestive system, including pancreatic secretion, as well as the digestion and absorption of vitamins and proteins in virulent BPHs (Figures 4B,D). These findings corroborate the importance of the digestive system’s importance in BPH virulence. In sum, the findings indicate that lipid metabolism and digestive system pathways are strongly involved in BPHs’ adaptation to resistant YHY15 plants.

Previous studies showed that the polygenic system controls BPH virulence (Claridge and Den Hollander, 1980; Tanaka, 1999; Jing et al., 2017). For example, BPH virulence to rice plants with Bph1-mediated resistance is influenced by a recessive gene located on linkage group 10 (Kobayashi et al., 2014), a gene on chromosome 7 influences the preference for Bph1 plants, and two major QTLs on chromosomes 5 and 14 influence the insect’s growth rates on resistant rice plants (Jing et al., 2014). The BSR-seq approach in this study extended these findings by identifying 24 SNPs linked with BPH virulence on six chromosomes. Most of these SNPs are located on chromosomes 3, followed by chromosomes 13, 4, 7, 8, and 14 (Supplementary Table 6). Furthermore, two chromosome regions (Chr3, 73413752-83702207 and 88874910-94058321) were identified (Figure 6). All the results suggest that the mechanism underlying BPH virulence is complex and involves multiple genes. In this analysis, we crossed a virulent female from biotype Y and an avirulent male from biotype 1 to produce an F2 population and performed a BSR-sequencing between the two extreme bulks. In the future, more BPH crosses and repeated BSR-sequencing analysis should be carried out to confirm and fine map the identified loci. Loci harboring the SNPs and candidate genes are expected to greatly facilitate determination of mechanisms involved in BPH virulence.

Conclusion

We identified 751 DEGs between the virulent and avirulent BPHs and mapped 24 SNPs linked to BPH virulence to rice plants with Bph15-mediated resistance by BSR-seq analysis. Functional analysis of 418 DEGs highly expressed in the avirulent insects showed that BPH upregulates carbohydrate, amino acid, and nucleotide metabolism, endocrine system, and signal transduction pathways in response to feeding on these plants. Lipid metabolism and digestive system pathways may play roles for BPH adaption to resistant rice. We expect the identification of relevant genetic loci and genes presented here to facilitate further research on BPH virulence mechanisms and responses to resistant plants.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA792102.

Author Contributions

WG and GH designed the research, analyzed data, and wrote the manuscript. WG performed most of the experiments. JS, MG, JG, DW, QZ, JW, RC, BD, LZ, and GH performed some of the experiments. All authors read and approved the manuscript.

Funding

This work was supported by a grant (no. 31630063) from the National Natural Science Foundation of China.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.843227/full#supplementary-material

Footnotes

1. ^ http://samtools.sourceforge.net/mpileup.shtml
2. ^ https://daehwankimlab.github.io/hisat2/

References

Abd-El-Haliem, A. M., Hoogstrate, S. W., and Schuurink, R. C. (2018). A robust functional genomics approach to identify effector genes required for thrips (Frankliniella occidentalis) reproductive performance on tomato leaf discs. Front. Plant Sci. 9:1852. doi: 10.3389/fpls.2018.01852

PubMed Abstract | CrossRef Full Text | Google Scholar

Abe, A., Kosugi, S., Yoshida, K., Natsume, S., Takagi, H., Kanzaki, H., et al. (2012). Genome sequencing reveals agronomically important loci in rice using MutMap. Nat. Biotechnol. 30, 174–178. doi: 10.1038/nbt.2095

PubMed Abstract | CrossRef Full Text | Google Scholar

Anders, S., and Huber, W. (2012). Differential Expression of RNA-Seq Data at the Gene Level – The DESeq Package. Heidelberg: European Molecular Biology Laboratory.

Google Scholar

Andersen, S. O. (2007). Involvement of tyrosine residues, N-terminal amino acids, and beta-alanine in insect cuticular sclerotization. Insect Biochem. Mol. Biol. 37, 969–974. doi: 10.1016/j.ibmb.2007.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Barski, O. A., Tipparaju, S. M., and Bhatnagar, A. (2009). Kinetics of nucleotide binding to the beta-subunit (AKR6A2) of the voltage-gated potassium (Kv) channel. Chem. Biol. Interact. 178, 165–170. doi: 10.1016/j.cbi.2008.10.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Bock, K. W. (2016). The UDP-glycosyltransferase (UGT) superfamily expressed in humans, insects and plants: animal-plant arms-race and co-evolution. Biochem. Pharmacol. 99, 11–17. doi: 10.1016/j.bcp.2015.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C. Y., Liu, Y. Q., Song, W. M., Chen, D. Y., Chen, F. Y., Chen, X. Y., et al. (2019). An effector from cotton bollworm oral secretion impairs host plant defense signaling. Proc. Natl. Acad. Sci. U.S.A. 116, 14331–14338. doi: 10.1073/pnas.1905471116

PubMed Abstract | CrossRef Full Text | Google Scholar

Chepelev, I., Wei, G., Tang, Q., and Zhao, K. (2009). Detection of single nucleotide variations in expressed exons of the human genome using RNA-Seq. Nucleic Acids Res. 37:e106. doi: 10.1093/nar/gkp507

PubMed Abstract | CrossRef Full Text | Google Scholar

Cingolani, P., Platts, A., Wang le, L., Coon, M., Nguyen, T., Wang, L., et al. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 6, 80–92. doi: 10.4161/fly.19695

PubMed Abstract | CrossRef Full Text | Google Scholar

Claridge, M. F., and Den Hollander, J. (1980). The “biotypes” of the rice brown planthopper, Nilaparvata lugens. Entomol. Exp. Appl. 27, 23–30. doi: 10.1111/j.1570-7458.1980.tb02942.x

CrossRef Full Text | Google Scholar

Claridge, M. F., and Den Hollander, J. (1982). Virulence to rice cultivars and selection for virulence in populations of the brown planthopper Nilaparvata lugens. Entomol. Exp. Appl. 32, 213–221. doi: 10.1111/j.1570-7458.1982.tb03208.x

CrossRef Full Text | Google Scholar

Danecek, P., Bonfield, J. K., Liddle, J., Marshall, J., Ohan, V., Pollard, M. O., et al. (2021). Twelve years of SAMtools and BCFtools. Gigascience 10:giab008. doi: 10.1093/gigascience/giab008

PubMed Abstract | CrossRef Full Text | Google Scholar

Den Hollander, J., and Pathak, P. K. (1981). The genetics of the “biotypes” of the rice brown planthopper, Nilaparvata lugens. Entomol. Exp. Appl. 29, 76–86. doi: 10.1111/j.1570-7458.1981.tb03044.x

CrossRef Full Text | Google Scholar

Du, H., Zhu, J., Su, H., Huang, M., Wang, H., Ding, S., et al. (2017). Bulked Segregant RNA-seq reveals differential expression and SNPs of candidate genes associated with waterlogging tolerance in maize. Front. Plant Sci. 8:1022. doi: 10.3389/fpls.2017.01022

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, X., Li, X., Yang, X., and Zhu, P. (2020). Fine mapping and identification of the leaf shape gene BoFL in ornamental kale. TAG 133, 1303–1312. doi: 10.1007/s00122-020-03551-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Foyer, C. H., Noctor, G., and van Emden, H. F. (2007). An evaluation of the costs of making specific secondary metabolites: does the yield penalty incurred by host plant resistance to insects result from competition for resources? Int. J. Pest Manage. 53, 175–182. doi: 10.1080/09670870701469146

CrossRef Full Text | Google Scholar

Górczak, K., Claesen, J., and Burzykowski, T. (2020). A conceptual framework for abundance estimation of genomic targets in the presence of ambiguous short sequencing reads. J. Comput. Biol. 27, 1232–1247. doi: 10.1089/cmb.2019.0272

PubMed Abstract | CrossRef Full Text | Google Scholar

Grbić, M., Van Leeuwen, T., Clark, R. M., Rombauts, S., Rouzé, P., Grbić, V., et al. (2011). The genome of Tetranychus urticae reveals herbivorous pest adaptations. Nature 479, 487–492. doi: 10.1038/nature10640

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J., Xu, C., Wu, D., Zhao, Y., Qiu, Y., Wang, X., et al. (2018). Bph6 encodes an exocyst-localized protein and confers broad resistance to planthoppers in rice. Nat. Genet. 50, 297–306. doi: 10.1038/s41588-018-0039-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hargrove, J. W. (1976). Amino acid metabolism during flight in tsetse flies. J. Insect Physiol. 22, 309–313. doi: 10.1016/0022-1910(76)90040-8

CrossRef Full Text | Google Scholar

He, P., Wei, P., Ma, Y., Hu, S., Yao, J., Jiang, X., et al. (2022). Candidate sex-associated gene identification in Trachinotus ovatus (Carangidae) using an integrated SLAF-seq and bulked segregant analysis approach. Gene 809:146026. doi: 10.1016/j.gene.2021.146026

PubMed Abstract | CrossRef Full Text | Google Scholar

Hibino, H. (1996). Biology and epidemiology of rice viruses. Annu. Rev. Phytopathol. 34, 249–274. doi: 10.1146/annurev.phyto.34.1.249

PubMed Abstract | CrossRef Full Text | Google Scholar

Hisano, H., Sakamoto, K., Takagi, H., Terauchi, R., and Sato, K. (2017). Exome QTL-seq maps monogenic locus and QTLs in barley. BMC Genomics 18:125. doi: 10.1186/s12864-017-3511-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Holtof, M., Lenaerts, C., Cullen, D., and Vanden Broeck, J. (2019). Extracellular nutrient digestion and absorption in the insect gut. Cell Tissue Res. 377, 397–414. doi: 10.1007/s00441-019-03031-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Jena, K. K., and Kim, S. M. (2010). Current status of brown planthopper (BPH) resistance and genetics. Rice 3, 161–171. doi: 10.1007/s12284-010-9050-y

CrossRef Full Text | Google Scholar

Ji, R., Yu, H., Fu, Q., Chen, H., Ye, W., Li, S., et al. (2013). Comparative transcriptome analysis of salivary glands of two populations of rice brown planthopper, Nilaparvata lugens, that differ in virulence. PLoS One 8:e79612. doi: 10.1371/journal.pone.0079612

PubMed Abstract | CrossRef Full Text | Google Scholar

Jing, S., Liu, B., Peng, L., Peng, X., Zhu, L., Fu, Q., et al. (2012). Development and use of EST-SSR markers for assessing genetic diversity in the brown planthopper (Nilaparvata lugens Stål). Bull. Entomol. Res. 102, 113–122. doi: 10.1017/S0007485311000435

PubMed Abstract | CrossRef Full Text | Google Scholar

Jing, S., Zhang, L., Ma, Y., Liu, B., Zhao, Y., Yu, H., et al. (2014). Genome-wide mapping of virulence in brown planthopper identifies loci that break down host plant resistance. PLoS One 9:e98911. doi: 10.1371/journal.pone.0098911

PubMed Abstract | CrossRef Full Text | Google Scholar

Jing, S., Zhao, Y., Du, B., Chen, R., Zhu, L., and He, G. (2017). Genomics of interaction between the brown planthopper and rice. Curr. Opin. Insect Sci. 19, 82–87. doi: 10.1016/j.cois.2017.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, T. (2016). Evolving ideas about genetics underlying insect virulence to plant resistance in rice-brown planthopper interactions. J. Insect Physiol. 84, 32–39. doi: 10.1016/j.jinsphys.2015.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, T., Yamamoto, K., Suetsugu, Y., Kuwazaki, S., Hattori, M., Jairin, J., et al. (2014). Genetic mapping of the rice resistance-breaking gene of the brown planthopper Nilaparvata lugens. Proc. Biol. Sci. 281:20140726. doi: 10.1098/rspb.2014.0726

PubMed Abstract | CrossRef Full Text | Google Scholar

Laity, J. H., Lee, B. M., and Wright, P. E. (2001). Zinc finger proteins: new insights into structural and functional diversity. Curr. Opin. Struct. Biol. 11, 39–46. doi: 10.1016/s0959-440x(00)00167-6

CrossRef Full Text | Google Scholar

Lam, S., Zeidan, J., Miglior, F., Suárez-Vega, A., Gómez-Redondo, I., Fonseca, P., et al. (2020). Development and comparison of RNA-sequencing pipelines for more accurate SNP identification: practical example of functional SNP detection associated with feed efficiency in Nellore beef cattle. BMC Genomics 21:703. doi: 10.1186/s12864-020-07107-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Li, D. L., Liu, S. Z., Ma, X. L., Dietrich, C. R., Hu, H. C., et al. (2013). The maize glossy13 gene, cloned via BSR-seq and seq-walking, encodes a putative ABC transporter required for the normal accumulation of epicuticular waxes. PLoS One 9:e99563. doi: 10.1371/journal.pone.0082333

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q., Pan, Z., Gao, Y., Li, T., Liang, J., Zhang, Z., et al. (2020). Quantitative Trait Locus (QTLs) mapping for quality traits of wheat based on high density genetic map combined with Bulked Segregant Analysis RNA-seq (BSR-Seq) indicates that the basic 7s globulin gene is related to falling number. Front. Plant Sci. 11:600788. doi: 10.3389/fpls.2020.600788

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C., Du, B., Hao, F., Lei, H., Wan, Q., He, G., et al. (2017). Dynamic metabolic responses of brown planthoppers towards susceptible and resistant rice plants. Plant Biotechnol. J. 15, 1346–1357. doi: 10.1111/pbi.12721

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S. Z., Yeh, C. T., Tang, H. M., Nettleton, D., and Schnable, P. S. (2012). Gene mapping via bulked segregant RNA-seq (BSR-seq). PLoS One 7:e36406. doi: 10.1371/journal.pone.0036406

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv, W., Du, B., Shangguan, X., Zhao, Y., Pan, Y., Zhu, L., et al. (2014). BAC and RNA sequencing reveal the brown planthopper resistance gene BPH15 in a recombination cold spot that mediates a unique defense mechanism. BMC Genomics 15:674. doi: 10.1186/1471-2164-15-674

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, W., Xu, L., Hua, H., Chen, M., Guo, M., He, K., et al. (2021). Chromosomal-level genomes of three rice planthoppers provide new insights into sex chromosome evolution. Mol. Ecol. Resour. 21, 226–237. doi: 10.1111/1755-0998.13242

PubMed Abstract | CrossRef Full Text | Google Scholar

Michelmore, R. W., Paran, I., and Kesseli, R. V. (1991). Identification of markers linked to disease-resistance genes by bulked segregant analysis: a rapid method to detect markers in specific genomic regions by using segregating populations. Proc. Natl. Acad. Sci. U.S.A. 88, 9828–9832. doi: 10.1073/pnas.88.21.9828

PubMed Abstract | CrossRef Full Text | Google Scholar

Navarro-Escalante, L., Zhao, C., Shukle, R., and Stuart, J. (2020). BSA-Seq discovery and functional analysis of candidate hessian fly (Mayetiola destructor) avirulence genes. Front. Plant Sci. 11:956. doi: 10.3389/fpls.2020.00956

PubMed Abstract | CrossRef Full Text | Google Scholar

Nelson, R. E., Fessler, L. I., Takagi, Y., Blumberg, B., Keene, D. R., Olson, P. F., et al. (1994). Peroxidasin: a novel enzyme-matrix protein of drosophila development. EMBO J. 13, 3438–3447. doi: 10.1002/j.1460-2075.1994.tb06649.x

CrossRef Full Text | Google Scholar

Novoselov, A., Becker, T., Pauls, G., von Reuß, S. H., and Boland, W. (2015). Spodoptera littoralis detoxifies neurotoxic 3-nitropropanoic acid by conjugation with amino acids. Insect Biochem. Mol. Biol. 63, 97–103. doi: 10.1016/j.ibmb.2015.05.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Otani, Y., Nakatsu, Y., Sakoda, H., Fukushima, T., Fujishiro, M., Kushiyama, A., et al. (2013). Integrator complex plays an essential role in adipose differentiation. Biochem. Biophys. Res. Commun. 434, 197–202. doi: 10.1016/j.bbrc.2013.03.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Pathak, P. K., and Heinrichs, E. A. (1982). Selection of biotype populations 2 and 3 of Nilaparvata lugens by exposure to resistant rice varieties. Environ. Entomol. 11, 85–90. doi: 10.1093/ee/11.1.85

CrossRef Full Text | Google Scholar

Pauchet, Y., Muck, A., Svatos, A., Heckel, D. G., and Preiss, S. (2008). Mapping the larval midgut lumen proteome of Helicoverpa armigera, a generalist herbivorous insect. J. Proteome Res. 7, 1629–1639. doi: 10.1021/pr7006208

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, L., Zhao, Y., Wang, H., Song, C., Shangguan, X., Ma, Y., et al. (2017). Functional study of cytochrome p450 enzymes from the brown planthopper (Nilaparvata lugens Stål) to analyze its adaptation to BPH-resistant rice. Front. Physiol. 8:972. doi: 10.3389/fphys.2017.00972

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez-Rubio, P., Lottaz, C., and Engelmann, J. C. (2019). FastqPuri: high-performance preprocessing of RNA-seq data. BMC Bioinformatics 20:226. doi: 10.1186/s12859-019-2799-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajan-Babu, I. S., Peng, J. J., Chiu, R., Imagine Study, Causes Study, Li, C., et al. (2021). Genome-wide sequencing as a first-tier screening test for short tandem repeat expansions. Genome Med. 13:126. doi: 10.1186/s13073-021-00932-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, A., Trapnell, C., Donaghey, J., Rinn, J. L., and Pachter, L. (2011). Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biol. 12:R22. doi: 10.1186/gb-2011-12-3-r22

PubMed Abstract | CrossRef Full Text | Google Scholar

Saxena, R. C., Demayo, C. G., and Barrion, A. A. (1991). Allozyme variation among biotypes of the brown planthopper Nilaparvata lugens in the Philippines. Biochem. Genet. 29, 115–123. doi: 10.1007/BF00554205

CrossRef Full Text | Google Scholar

Shi, S., Wang, H., Nie, L., Tan, D., Zhou, C., Zhang, Q., et al. (2021). Bph30 confers resistance to brown planthopper by fortifying sclerenchyma in rice leaf sheaths. Mol. Plant. 14, 1714–1732. doi: 10.1016/j.molp.2021.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Smadja, C., Shi, P., Butlin, R. K., and Robertson, H. M. (2009). Large gene family expansions and adaptive evolution for odorant and gustatory receptors in the pea aphid, Acyrthosiphon pisum. Mol. Biol. Evol. 26, 2073–2086. doi: 10.1093/molbev/msp116

PubMed Abstract | CrossRef Full Text | Google Scholar

Sogawa, K. (1981). Hybridization experiments on 3 biotypes of the brown planthopper, Nilaparvata lugens (Homoptera: Delphacidae) at the IRRI, the Philippines. Appl. Entomol. Zool. 16, 193–199. doi: 10.1303/aez.16.193

PubMed Abstract | CrossRef Full Text | Google Scholar

Sogawa, K. (1982). The rice brown planthopper: feeding physiology and host plant interactions. Ann. Rev. Entomol. 27, 49–73. doi: 10.1146/annurev.en.27.010182.000405

CrossRef Full Text | Google Scholar

Srinivasan, K. A., Virdee, S. K., and McArthur, A. G. (2020). Strandedness during cDNA synthesis, the stranded parameter in htseq-count and analysis of RNA-Seq data. Brief. Funct. Genomics 19, 339–342. doi: 10.1093/bfgp/elaa010

PubMed Abstract | CrossRef Full Text | Google Scholar

Takagi, H., Abe, A., Yoshida, K., Kosugi, S., Natsume, S., Mitsuoka, C., et al. (2013). QTL-seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations. Plant J. 74, 174–183. doi: 10.1111/tpj.12105

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanaka, K. (1999). Quantitative genetic analysis of biotypes of the brown planthopper Nilaparvata lugens: heritability of virulence to resistant rice varieties. Entomol. Exp. Appl. 90, 279–287. doi: 10.1046/j.1570-7458.1999.00448.x

CrossRef Full Text | Google Scholar

Teerawanichpan, P., Robertson, A. J., and Qiu, X. (2010). A fatty acyl-CoA reductase highly expressed in the head of honey bee (Apis mellifera) involves biosynthesis of a wide range of aliphatic fatty alcohols. Insect Biochem. Mol. Biol. 40, 641–649. doi: 10.1016/j.ibmb.2010.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Tereshchenkova, V. F., Goptar, I. A., Kulemzina, I. A., Zhuzhikov, D. P., Serebryakova, M. V., Belozersky, M. A., et al. (2016). Dipeptidyl peptidase 4 – An important digestive peptidase in Tenebrio molitor larvae. Insect Biochem. Mol. Biol. 76, 38–48. doi: 10.1016/j.ibmb.2016.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621

PubMed Abstract | CrossRef Full Text | Google Scholar

Via, S. (1990). Ecological genetics and host adaptation in herbivorous insects: the experimental study of evolution in natural and agricultural systems. Annu. Rev. Entomol. 35, 421–446. doi: 10.1146/annurev.en.35.010190.002225

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, B., Shahzad, M. F., Zhang, Z., Sun, H., Han, P., Li, F., et al. (2014). Genome-wide analysis reveals the expansion of Cytochrome P450 genes associated with xenobiotic metabolism in rice striped stem borer, Chilo suppressalis. Biochem. Biophys. Res. Commun. 443, 756–760. doi: 10.1016/j.bbrc.2013.12.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R. J., Sun, L. Y., Bao, L. S., Zhang, J. R., Jiang, Y. L., Yao, J., et al. (2013). Bulk segregant RNA-seq reveals expression and positional candidate genes and allele-specific expression for disease resistance against enteric septicemia of catfish. BMC Genomics 14:929. doi: 10.1186/1471-2164-14-929

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Zhang, M., Feng, F., and He, R. (2015). Differentially regulated genes in the salivary glands of brown planthopper after feeding in resistant versus susceptible rice varieties. Arch. Insect Biochem. Physiol. 89, 69–86. doi: 10.1002/arch.21226

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, H., You, A., Yang, Z., Zhang, F., He, R., Zhu, L., et al. (2004). High-resolution genetic mapping at the Bph15 locus for brown planthopper resistance in rice (Oryza sativa L.). Theor. Appl. Genet. 110, 182–191. doi: 10.1007/s00122-004-1844-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Yue, L., Kang, K., and Zhang, W. (2019). Metabolic responses of brown planthoppers to IR56 resistant rice cultivar containing multiple resistance genes. J. Insect Physiol. 113, 67–76. doi: 10.1016/j.jinsphys.2018.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Guan, W., Huang, C., Hu, Y., Chen, Y., Guo, J., et al. (2019). Combining next-generation sequencing and single-molecule sequencing to explore brown plant hopper responses to contrasting genotypes of japonica rice. BMC Genomics 20:682. doi: 10.1186/s12864-019-6049-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, X., Xin, Y., Peng, Y., Shan, J., Zhang, N., Wu, D., et al. (2020). Lipidomic analyses reveal enhanced lipolysis in planthoppers feeding on resistant host plants. Sci. China Life Sci. 64, 1502–1521. doi: 10.1007/s11427-020-1834-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, X., Zhu, L., and He, G. (2021). Genetic and molecular understanding of host rice resistance and Nilaparvata lugens adaptation. Curr. Opin. Insect Sci. 45, 14–20. doi: 10.1016/j.cois.2020.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, C., Zhang, Q., Chen, Y., Huang, J., Guo, Q., Li, Y., et al. (2021). Balancing selection and wild gene pool contribute to resistance in global rice germplasm against planthopper. J. Integr. Plant Biol. 63, 1695–1711. doi: 10.1111/jipb.13157

PubMed Abstract | CrossRef Full Text | Google Scholar

Zielinski, T., Reichman, M., Donover, P. S., and Lowery, R. G. (2016). Development and validation of a universal high-throughput UDP-glycosyltransferase assay with a time-resolved FRET signal. Assay Drug Dev. Technol. 14, 240–251. doi: 10.1089/adt.2016.711

PubMed Abstract | CrossRef Full Text | Google Scholar

Zogli, P., Pingault, L., Grover, S., and Louis, J. (2020). Ento(o)mics: the intersection of ‘omic’ approaches to decipher plant defense against sap-sucking insect pests. Curr. Opin. Plant Biol. 56, 153–161. doi: 10.1016/j.pbi.2020.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar