Identification of the Carboxylesterase, Nl-EST1 Gene Expression Regulation Possibly Involved in Carbamate Insecticide Resistance in Brown Planthopper Nilaparvata lugens

1. Introduction

Rice is widely consumed staple grain rich in fiber, energy, minerals, vitamins, and various biomolecules [1]. Research has shown that different parts of rice offer numerous health benefits in both preclinical and clinical studies. As a result, the constituents of rice are becoming increasingly popular for use in creating pharmaceutical adjuvants, food additives, and dietary supplements [2]. However, rice is a suitable host for a wide range of insects. Insect pests have always been a major issue for rice, causing significant yield loss and deterioration in grain quality. Among these pests, the brown planthopper (BPH, Nilaparvata lugens Stål, 1854) (Hemiptera: Delphacidae) is the most destructive, leading to about 20% to 80% yield loss and an annual economic loss of around $300 million in Asia [3]. BPH harms rice crops by extracting sap from the xylem and phloem tissues, resulting in ‘hopper burn’. In addition, BPH indirectly causes harm by transmitting viral diseases like grassy stunt virus and ragged stunt virus [4,5].

In Korea, N. lugens, a significant migratory insect species originating from China, seriously threatens agriculture [6]. This pest is primarily controlled using organophosphate (OP), carbamate (CB), pyrethroid, nereistoxin and neonicotinoid insecticides. OP and CB insecticides have been employed for 50 years to manage N. lugens populations. However, the extensive use of OP and CB insecticides has led N. lugens to develop resistance in Japan, Taiwan, Solomon Islands, Philippines, Malaysia and Korea [7].

Metabolic resistance in insects to insecticides involves the detoxification of chemicals through the over-production of specific enzymes that break down the insecticides before they can reach and bind to their target sites. This mechanism relies heavily on enzymes such as monooxygenases (cytochrome P450 monooxygenases, CYP), hydrolases (such as esterases, ESTs), and transferases (glutathione-S-transferase, GST), which convert xenobiotics into non-toxic compounds [8]. Esterases, such as E4 and FE4, degrade ester bonds in insecticides like OP and CB, often due to gene amplification or upregulation [9]. Overexpression of CYPs, particularly in species like Myzus persicae is linked to neonicotinoid resistance through increased gene copy numbers and mutations [10]. Similarly, GSTs are pivotal in detoxifying endogenous and exogenous compounds, with enhanced expression linked to resistance to various pests [11]. Mutation-based resistance occurs when the target insecticide site is modified, reducing binding efficiency. Notable examples include nicotinic acetylcholine receptors where single point mutations in the D-loop region or alterations in subunits like Mdα2 and Mdα6 in houseflies confer resistance to neonicotinoids and spinosad [12,13,14]. Additionally, modified acetylcholine esterase due to gene mutations results in insensitivity to OP and CB, with specific mutations like G262V in Musca domestica showing strong resistance [15]. Both metabolic and mutation-based resistance exemplify pests' adaptive strategies to overcome chemical control measures, necessitating ongoing research and development of novel insecticides.

Carboxyl/choline esterase (CCE), particularly carboxylesterases constitute a diverse and widespread group of enzymes involved in a range of metabolic processes, including hormone metabolism, pheromone breakdown, detoxification of foreign substances, and hydrolysis of carboxyl esters in insecticides [16]. The enhanced ability of CCEs to detoxify is linked to resistance to various insecticides like organophosphates, carbamates, and pyrethroids [17]. Previous research indicates that heightened expression and activity of specific CCE are associated with increased resistance to insecticides across several insect species, such as N. lugens [18] and Myzus persicae [19,20].

Our research group identified the I4790M mutation in the Spodoptera exigua ryanodine receptor gene, noting that its correlation with resistance to diamide insecticides varied among the field-collected population in Korea [21]. The authors proposed that transcription regulation could contribute to enhanced resistance in S. exigua against diamides in addition to this mutation. Similarly, Adelman, et al. [22] identified that the esterase-encoding genes CE3959 and CE21331 were significantly overexpressed in the highly resistant Richmond strain of Cimex lectularius, suggesting their role in esterase-mediated resistance. Zhu, et al. [23] also confirmed the overexpression of CLCE21331 in resistance strains of C. lectularius, showing more than 50-fold up-regulation in most field populations, indicating its critical role in pyrethroid resistance. In contrast, mutation in coding gene sequences is rarely reported globally and has not yet been observed in bed bugs, indicating its limited documentation and potential relevance [24]. Building on this, our current study aimed to investigate resistance to the carbamate insecticide fenobucarb in N. lugens in Korea, focusing on transcriptomic changes.

2. Materials and Methods

2.1. Insects

The insecticide-susceptible strain (BPH80) of brown planthopper (BPH, Nilaparvata lugens) was shared by the National Academy of Agriculture Science, Rural Development Administration (RDA), Korea. This strain was collected from the field in 1980 and grown in the laboratory without insecticide exposure for over 40 years. In addition to the insecticide-susceptible strain collected in 1980, they were collected from Wando (34°22'06"N, 126°43'27"E) and Goseong (34°57'19"N, 128°20'38"E) in 2015 (BPH15) and 2019 (BPH19), respectively, and reared in the National Institute of Crop Science, RDA. All bioassays and RNA extractions were performed after 2019 when BPH19 was secured, there was no more selection by insecticide treatment in the lab. BPH80 was maintained for about 40 years, and BPH15 for about five years without insecticide exposure in the lab. All experiments were performed on BPH19 within three generations after collection. All BPHs were maintained under controlled conditions at a temperature of 25 ± 1 °C, relative humidity of 60 ± 5%, and a photoperiod regime of 14 hours of light and 10 hours of darkness as previously reported [7].

2.2. Bioassay

Bioassays were performed based on the IRAC (Insecticide Resistance Action Committee) susceptibility test method 005 with some modifications (www.irac-online.org) for adults. Fenobucarb (BPMC, Emulsifiable Concentrate Formulation of 50 % w/v) was used for bioassay. After diluting fenobucarb to various concentrations, rice seedlings ten days after seeding were dipped into the diluted solution for 10 seconds, and then ten adults per treatment were transferred. All experiments conducted more than three biological replications (n>30 per concentration). More detailed bioassay methods were followed by IRAC susceptibility test method 005.

Using the SAS program based on the Probit model (SAS Institute 9.1, Cary, NC, USA), concentration-based mortality after two days of fenobucarb exposure was estimated to determine the median lethal concentration (LC₅₀) and 95% confidence limits (CLs). RR (resistance ratio) was computed by dividing the LC₅₀ value of the tested field population by that of the susceptible strain, BPH80.

2.3. RNA and DNA Extraction

Total RNAs were extracted from the adults of each strain and populations of N. lugens within 12 h after emergence, with each sample containing twenty adults as a biological replicate. RNA extraction was conducted using the RNeasy Mini Kit (Qiagen, Hilden, Germany), according to the manufacturer’s instructions. The RNA was validated and quantified using an Agilent 2200 TapeStation (Agilent Technologies, Santa Clara, CA, USA), and RNA integrity was confirmed by running samples on a 1% agarose gel using electrophoresis. For the reverse transcription reaction, we utilized the SuperiorScript III cDNA Synthesis Kit (Enzynomics, Daejeon, Korea). The total RNAs and synthesized cDNA were stored at −70 °C before the next experiments. Furthermore, genomic DNA (gDNA) was extracted from the adults of each N. lugens within 12 h after emergence, with each sample containing twenty adults as a biological replicate using DNeasy Blood & Tissue (Qiagen) following the manufacturer’s instructions and quantified using Nanodrop (Nanodrop Technologies, Wilmington, DE, USA).

2.4. Mutation Survey

Following specific thermal conditions, the cDNA and gDNA underwent PCR to survey the mutations in ace1 and Nl-EST1, using the ProFlex PCR System (ThermoFisher Scientific) with KOD FX polymerase (Toyobo Life Science, Osaka, Japan) with the appropriate primer combinations and PCR conditions. Used primer sets are listed in Table 1. The PCR products were directly sequenced (Macrogen, Seoul, Korea) and the chromatograms were analyzed for mutations, following the methodology described earlier [21].

2.5. RNA-seq Analysis

Total RNAs were extracted from the adults of a strain (BPH80) and populations (BPH15 and BPH19) of N. lugens within 12 h after emergence, with each sample containing twenty adults as a biological replicate. RNA extraction was conducted using the RNeasy Mini Kit (Qiagen, Hilden, Germany), according to the manufacturer’s instructions. The RNA was validated and quantified using an Agilent 2200 TapeStation (Agilent Technologies, Santa Clara, CA, USA), and RNA integrity was confirmed by running samples on a 1% agarose gel using electrophoresis. RNA-seq libraries were prepared using the TruSeq RNA sample Prep Kit v2 (Illumina, San Diego, CA, USA). Samples were sequenced on the Hiseq4000 plDEGatform using TruSeq 3000/4000 SBS Kit v3 (Macrogen, Seoul, Korea).

The nine RNA-seq raw sequences were initially processed using Trimmomatic v0.38 to remove low-quality sequences (Q30) and adapters from the raw data [25]. The quality of the resulting trimmed reads was confirmed to be high using FastQC v0.11.7 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc).

2.6. Clean Read Assembly and Unigene Construction

Trimmed reads from all samples were merged into each assembly group to construct transcriptome references. The merged data was assembled using the Trinity version (r20140717) program, utilized for de novo transcriptome assembly [26]. This process results in transcript fragments called contigs. The longest contigs were clustered into non-redundant transcripts, referred to as unigenes, using the CD-HIT-EST program provided by CD-HIT v4.6 [27].

2.7. Functional Annotation

Protein sequences were generated by predicting ORFs from unigenes using TransDecoder v3.0.1, and the generated unigenes and protein sequences were used for annotation and expression analysis. To perform functional annotation, we used BLASTN from NCBI BLAST v2.9.0+ and BLASTX of DIAMOND v0.9.21 software with an E-value default cutoff of 1.0E-5 [28], against the Kyoto Encyclopedia of Genes and Genomes (KEGG), NCBI Nucleotide (NT), Pfam, Gene Ontology (GO), NCBI non-redundant Protein (NR), UniProt, and EggNOG [29].

2.8. Differential Gene Expression (DGE) Analysis

Trimmed reads for each sample were aligned to the assembled unigene as reference using Bowtie v1.1.2. For the differentially expressed gene analysis, the abundances of unigenes across samples are estimated into the read count as an expression measure by the RSEM v1.3.1 algorithm [30]. For nine samples, if more than one read count value was 0, it was not included in the analysis. To reduce systematic bias, estimate the size factors from the count data and apply Relative Log Expression (RLE) normalization with DESeq2 v1.28.1 (https://www.bioconductor.org/packages/release/bioc/html/DESeq2.html). Differentially expressed genes (DEGs) analysis and statistical analysis were performed using Log2FoldChange (FC), nbinomWaldTest per comparison pair using DESeq2 (|FC| >2 and nbinomWaldTest raw p-value<0.05).

2.9. Orthologous Cluster Analysis

Using the web program Orthovenn3, orthologous and particular genes among the predicted protein unigenes of all BPH were investigated and displayed [31]. Venn diagrams were used for the comparative analysis and visualization of orthologous gene clusters and distinct genes that are particular to each BPH. Then, protein sequences classified into strain-specific clusters were subjected to GO analysis to evaluate strain-specific biological roles.

3. Results

3.1. Bioassay

The rice seedling dip method-based bioassays demonstrated significant variations in susceptibility to the carbamate insecticide fenobucarb among the susceptible strain (BPH80), and field populations (BPH15 and BPH19). The LC₅₀ values were 3.08 for the BPH80, 10.61 for the BPH15 population, and 73.98 for the BPH19 population. Compared to BPH15, and BPH19 based on BPH80, exhibited resistance levels 3.4 times and 24.2 times higher, respectively. This indicates a progressive increase in resistance, with the order of susceptibility being BPH80>BPH15>BPH19.

These results highlight the growing resistance to fenobucarb and suggest that the BPH80 is highly susceptible, while BPH15 and BPH19 populations have developed significant resistance.

3.2. Mutation Analysis

The comparison of amino acid sequence in the mutation-reported regions of the ace1 gene across four pest species, including N. lugens, Laodelphax striatella, Sogatella furcifera (Delphacidae), and Nephotettix cincticeps (Cicadellidae), showed a highly conserved region with minimal variation. The ace2 gene in N. lugens was used as an outgroup for the analysis (Figure 1A). As a result, no mutations were found in the known resistance-associated sites within the field populations (Figure 1B, and C). These findings suggest that the fenobucarb resistance in these populations is not due to point mutations in the ace1 gene but may be related to transcriptomic changes of other genes. We checked for mutations in the Nl-EST1, but no nonsynonymous substitutions were found in BPH15 and BPH19 (Figure S1).

3.3. Raw and Trimmed Data Statistics of RNAseq

To ensure the reliability and quality of sequencing data for genetic analysis, we performed sequencing on three distinct samples: BPH80, BPH15, and BPH19, with three biological replicates each. The quality of raw sequencing data was evaluated by determining the total number of bases, total reads, GC contents (%), and Q30 (%) for nine samples. The total read bases per sample ranged from 5.39 giga base pair (Gb) to 6.46 Gb and total reads per sample ranged from 53.38 million to 64.03 million. The GC (%) content per sample ranged from 31 to 40. Additionally, Q30 (%) values were 95 and 96 as illustrated in Table S1. The data indicates high-quality sequencing across all samples, with consistently high Q30 values, ensuring reliable downstream analysis.

Subsequently, the raw reads were processed using the Trimmomatic program to remove adapter sequences and low-quality bases. The total read bases per sample ranged from 5.26 Gb to 6.33 Gb, and the total reads per sample ranged from 52.42 million to 63.10 million (Table S2). The GC (%) content per sample ranged from 31 to 40. Similarly, Q30 (%) values were 96 as illustrated in Table S2. This high-quality trimmed data ensures accuracy and reliability for subsequent analysis.

3.4. De Novo Assembly of Unigene Sets

To accurately represent gene expression patterns and reduce any potential errors, we carried out a de novo assembly analysis on samples from BPH80, BPH15, and BPH19 following preprocessing. The purpose of this analysis was to create comprehensive transcriptomes for each sample. Detailed statistics for the initial assembled contig and unigene contig are provided in Table 2 and Table S3, presenting data such as the number of genes, number of transcripts, %GC content, N50, average contig length, and total assembled bases (bp).

For the initial assembly, the emerged dataset revealed the highest number of genes at 148,234. In comparison, BPH19 had 111,118 genes, BPH80 had 92,362, and BPH15 had 79,691. The merged dataset also showed the highest number of transcripts at 191,287, surpassing BPH19 with 139,372, BPH80 with 115,494, and BPH15 with 97,626. The %GC content was highest in both the merged dataset and BPH19, both at 41%, while BPH80 and BPH15 had 39% and 38%, respectively. The N50 values mirrored this trend, with the merged assembly achieving an N50 value of 912, followed by BPH19 at 893, BPH80 at 810, and BPH15 with a lower value. The average contig length (bp) was greatest in the merged dataset at 626 bp, with BPH19 at 618 bp, BPH80 at 584 bp, and BPH15 at 577 bp. The total assembled bases were also highest in the merged dataset with 119,808,296 bp compared to 86,124,462 bp for BPH19, 67,502,190 bp for BPH80, and 56,350,010 bp for BPH15 (Table S3).

In the unigene contig assembly, the emerged dataset contained 119,664 genes, which was more than BPH19 with 93,427 genes, BPH80 with 75,069, and BPH15 with 69,319. The number of transcripts followed a similar pattern, with the merged dataset having 119,664 transcripts, BPH19 with 93,427, BPH80 with 75,069, and BPH15 with 69,319. The %GC content remained highest in both the merged dataset and BPH19 at 41%, with BPH80 at 39% and BPH15 at 38%. The N50 values for unigene contig were also highest in the merged assembly at 896, with BPH19 at 860, BPH80 at 816, and BPH15 with the lowest value. The average contig length (bp) was greatest in the merged assembly at 611 bp, followed by BPH19 at 600 bp, BPH80 at 579 bp, and BPH15 at 560 bp. Regarding total assembled bases, the merged dataset had 73,153,339 bp, while BPH19 had 56,127,492 bp, BPH80 had 43,504,047 bp, and BPH15 had 38,844,156 bp (Table S4).

3.5. ORF prediction

The analysis conducted on unigenes assembled from various datasets including merge, BPH80, BPH15, and BPH19, revealed notable variations in the proportion of unigenes containing predicted ORFs across these datasets. In the merged dataset, comparing 119,664 total unigenes, 25.73% (30,788 unigenes) were found to have predicted ORFs. It is worth mentioning that the majority of these, 97.65% (30,064 unigenes) were predicted to have a single ORF, while only 2.35% (724 unigenes) had multiple ORFs (Table 3 and Table S4).

In the BPH80 dataset, which consists of 75,069 total unigenes, 24.77% (18,593 unigenes) were identified as having predicted ORFs. Among these, 98.69% (18,349 unigenes) contained a single ORF, whereas a smaller fraction, 1.31% (244 unigenes), exhibited multiple ORFs. Similarly, the BPH15 dataset, with 69,319 total unigenes, showed that 22.91% (15,880 unigenes) were found to have predicted ORFs, of which 99.0% (15,721 unigenes) contained a single ORF, and only 1.0% (159 unigenes) had multiple ORFs. The BPH19 dataset, comprising 93,427 unigenes, demonstrated the highest percentage of predicted ORFs at 27.99% (26,147 unigenes). Among these, 98.7% (25,721 unigenes) had a single ORF, while 1.63% (426 unigenes) contained multiple ORFs (Table 3).

3.6. Transcriptomic DEG analysis

RNA-seq data from all samples were analyzed based on the merged reference unigene set as well as strain and population-specific unigenes. The mapping efficiency of reads to the reference unigene set (merged) was evaluated for various samples, providing insight into the alignment quality and comprehensiveness of the unigene assembly (Table S4). These results reflect the effectiveness of the unigene assembly in capturing the majority of the sequence reads, although a significant proportion of reads in some samples remained unmapped (Table S4).

The results of our hierarchical clustering unveiled varying susceptibility among the field populations (BPH15 and BPH19), and BPH80. BPH80 displayed higher susceptibility than BPH15 and BPH19 populations, which is illustrated in Figure 2A. Moreover, a principal component analysis (PCA) plot was created based on read counts from aligning each strain and the population’s sample reads to the reference unigene. The PCA analysis indicates distinct gene expression patterns for all BPHs (Figure 2B). Similarly, volcano plots revealed a significant increase in gene expression levels in both BPH15 and BPH19 compared to BPH80, as shown in Figure 2C and D. These RNA-seq analyses suggest that the upregulation of genes present in BPH15 and BPH19 populations conferred fenobucarb resistance compared to BPH80.

3.7. Orthologous Cluster Analysis

To gain a deeper understanding of the factors behind fenobucarb resistance, we analyzed the cluster count, protein count, singletons, and the number of specific or shared elements among all BPH. This analysis aimed to uncover the variations in gene expression and genetic diversity that may influence different levels of resistance. Following exposure to fenobucarb, we noted significant disparities in cluster counts. The BPH19 population, which exhibited higher resistance, showed a notably higher cluster count (16,708) compared to moderately resistant BPH15 (15,059) and the susceptible BPH80 (12,593), as illustrated in Figure 3A. Similarly, BPH19 demonstrated a higher number of clusters, proteins, and singletons compared to BPH15 and BPH80. Specifically, BPH19 had 16,708 clusters, 26,591 proteins, and 9,141 singletons. In contrast, BPH15 had 12,593 clusters, 16,045 proteins, and 3,030 singletons, while BPH80 had 15,059 clusters, 18,842 proteins, and 3,398 singletons, as shown in Figure 3B. Furthermore, we examined the presence of specific and shared elements among all BPHs. BPH19 displayed 454 specific elements, significantly more than BPH15 (21 specific elements) and BPH80 (20 specific elements). The analysis of shared elements revealed that BPH19 and BPH15 shared a total of 11645 elements (2142+9503), BPH19 and BPH80 shared a total of 14,112 elements (9,503+4,609), and BPH15 and BPH80 shared 10,430 elements (9503+927), as depicted in Figure 3C. By comparing genetic homology based on orthologous genes, BPH19, which had the highest resistance, was more similar to BPH80 than BPH15 (Figure 3D and E). Additional GO analysis was performed on 20 genes of BPH80, 21 genes of BPH15, and 454 genes of BPH19 expressed in a strain or population-specific manner in Figure 3C. These genes were categorized into three gene ontology (GO) categories: biological process, cellular component, and molecular function, as shown in Figure S2. Our analysis reveals that, in comparison to BPH80 and BPH15, BPH19 exhibited a higher representation across all three GO categories. Specifically, BPH19 demonstrated an enhanced number of unigenes associated with biological processes, cellular components and molecular functions, indicating a broader functional diversity in this genotype. The correlation matrix data revealed a strong correlation among all replicates of each BPH, which further confirms the accuracy and reliability of our results (Figure S3). This consistency across replicates underscores the validity of our findings related to gene expression and insecticide resistance.

3.8. DEG Analysis of Detoxification Genes

Differential gene expression (DGE) analysis was focused on the represented five detoxification gene families such as CYP, CCE, GST, uridine 5’-diphosphate-glucosyltransferase (UGT), ATP-binding cassette transporter (ABC), and cuticular protein (CP). Specifically, no correlation was found between resistance levels and the expression of CYPs, GSTs, UGTs, and ABCs in all strain and populations (Figure 4). Conversely, a strong correlation was observed between resistance levels and the expression of CCEs in PBH15, especially in BPH19, compared to BPH80 (Figure 4B).

To further confirm the correlation between the expression of CCEs and the resistance levels, we calculated fragment per kilobase of transcript per million map reads (FPKM) values using the number of RNA-Seq reads mapped to gene sequences for gene expression profiling. We examined the FPKM values to assess the expression of several genes from the CCE6 subfamily, E4 Type CCE, and CCE2 subfamily. The primary objective of this analysis was to identify the key gene or genes involved in enhancing resistance in N. lugens against fenobucarb (Figure 5). Our results revealed a highly significant increase in the expression of E4 type CCE gene NI-EST1 in BPH19, followed by BPH15, compared to BPH80 after fenobucarb exposure, as illustrated in Figure 5A. Conversely, no significant changes were observed in the expression of the CCE6 and CCE2 subfamily genes among all BPHs (Figure 5A). These findings indicate a strong correlation between the expression of the E4 type CCE gene NI-EST1 and heightened resistance levels. Resistance in BPHs has significantly increased, with BPH19 exhibiting the highest resistance to fenobucarb, followed by BPH15, while BPH80 remains highly sensitive to fenobucarb resistance.

4. Discussion

Rice is an essential staple food consumed worldwide, especially in Asia. It is also a common host for various pests, such as the brown planthopper N. lugens. This insect is highly destructive to rice plants as it feeds on their sap, leading to significant damage and the transmission of viral diseases, resulting in reduced crop yields. Plants have developed several defense mechanisms, including antibiosis, antixenosis, and tolerance to combat insect damage [32]. Additionally, the use of different insecticides can protect plants from the harmful effects of various insects [33]. However, the widespread and repeated use of insecticides has contributed to a notable increase in insect resistance, posing a significant challenge for formers and scientists globally [34].

Hence, we examined the rise in resistance of N. lugens collected in Korea in 1980, 2015 and 2019 to the carbamate insecticide fenobucarb. According to our bioassay results, BPH80 displayed greater susceptibility to fenobucarb than BPH15 and BPH19. Particularly, BPH19 exhibited markedly higher resistance to fenobucarb than both BPH80 and BPH15. The LC50 values were 3.08 for the BPH1980, 10.61 for the BPH2015, and 73.98 for the BPH2019. In comparison to the susceptible strain, the 2015 population demonstrated a resistance level of 3.4 times higher, while the 2019 population showed a resistance level of 24.2 times higher. These results align with prior studies that have noted an escalation in N. lugens resistance to carbamate [7,35,36]. The difference in the resistance level between BPH15 and BPH19 may be due to regional differences, but in the case of BPH15, it is also possible that the expression of resistance-related genes was downregulated due to prolonged lab rearing.

It appears that the ace1 gene conserved regions do not show any mutations across various pest species, indicating that fenobucarb resistance may not be caused by genetic changes at these loci in tested N. lugens samples (Figure 1). This discovery is consistent with previous research that has also noted the absence of mutation-related resistance, suggesting that alternative resistance mechanisms, such as transcriptomic changes, may play a critical role. These transcriptomic changes have been identified in other studies as important factors in modifying gene expression and regulatory pathways [21]. Further exploration of these transcriptomic mechanisms is essential to gain a comprehensive understanding of the observed resistance in these populations.

In contrast, specific point mutations were found in the ace1 gene of N. lugens, particularly F331H and I332L, that were associated with increased resistance levels to carbamate insecticides [35]. Their research indicated that these mutations could be used as molecular markers to detect resistance. However, our study did not find similar mutations in the conserved regions of ace1 across multiple pest species, including N. lugens. This disparity suggests that the resistance mechanisms in the populations we examined may be different, potentially involving alternative transcriptomic changes rather than point mutations. These differing outcomes highlight the complexity of insecticide resistance and emphasize the necessity for further research to investigate non-genetic factors that contribute to resistance in these pest populations.

To explore the reasons for the increased resistance to fenobucarb in BPH15 and BPH19 compared to BPH80, we utilized RNA-seq analysis. Our findings from hierarchical clustering (Figure 2A), multidimensional scaling (Figure 2B), and volcano plots (Figure 2C and D) indicate that the increased resistance observed in BPH15 and BPH19 is associated with significant changes in transcript accumulation. These results are consistent with previous studies suggesting that increased insecticide resistance in N. lugens is often due to elevated levels of transcript accumulation in resistant populations compared to susceptible ones [37]. Additionally, Pu, et al. [38] reported that multiple cis-acting elements contribute to the up-regulation of CYP6FU1, leading to resistance in L. striatellus. The higher resistance level of BPH19 compared to BPH15 and BPH80 was further confirmed on a genetic basis using the orthovenn 3 program. Our results indicated that BPH19 exhibits higher cluster counts, more clusters, proteins, and singletons than BPH15 and BPH80. These findings suggest that the increased resistance observed in BPH19 could be attributed to its greater genetic diversity compared to BPH15 and BPH80 (Figure 3A-C). Our results further support the notion that genetic diversity plays a crucial role in resistance mechanisms, highlighting how environmental stressors can drive phenotypic changes that enhance survival in pest populations. Malathi, et al. [39] also reported on the role of detoxifying enzymes in field-evolved resistance to various insecticides in N. lugens, underscoring the significance of metabolic resistance mechanisms. Notably, BPH19 exhibited a higher enrichment across all GO categories such as biological process, cellular component, and molecular function compared to BPH80 and BPH15, suggesting a greater functional diversity in BPH19 (Figure S2). Together these studies support the idea that genetic diversity and the upregulation of detoxifying enzymes are key factors in the development of insecticide resistance, reinforcing the robustness of our RNA-seq analysis in identifying genes in fenobucarb resistance.

The detoxification enzymes, such as CYPs, ESTs, GSTs, UGTs, ABCs, and CPs, are crucial in boosting insecticide resistance [40]. These enzymes have differential roles in regulating insect resistance to various insecticides. Our study indicates that when exposed to fenobucarb, the expression levels of the EST and CP gene families were significantly higher in BPH19 compared to BPH15 and BPH80, potentially contributing to the increased resistance of BPH19 (Figure 4). Rather than judging that these results contradict previous reports on mutation-based resistance development it was suggested that metabolism-based resistance development is also possible [35]. Different resistance development mechanisms can elevate the resistance of N. lugens when exposed to various insecticides [35]. Previously it was reported that the E4 and FE4 esterase genes enhance resistance in M. persicae within the esterase gene family [41]. Consequently, we also analyzed the FPKM values to evaluate the expression of various genes from the CCE6 subfamily, E4 Type CCE, and CCE2 subfamily etc. Our analysis uncovered a strong correlation between the resistance level of BPH19 and the expression of the EST1 gene, compared to BPH15 and BPH80 (Figure 5). Overall, our study suggests that the Nl-EST1 is strongly correlated with resistance of BPH to fenobucarb in N. lugens and may be a key gene in enhancing this resistance. Nonetheless, further research is necessary to validate these results and better understand the underlying mechanisms.

5. Conclusions

The use of carbamate insecticides such as fenobucarb to manage the rice pest BPH has been ongoing since their registration in Korea in 1981. Unfortunately, the improper and excessive application of these insecticides has resulted in heightened resistance in N. lugens populations, presenting a considerable challenge for farmers and researchers seeking to comprehend and address this resistance. Our bioassay analysis indicated that BPH19 exhibited significantly greater resistance than BPH80 and BPH15. To delve into the reasons behind BPH19’s enhanced resistance to fenobucarb, we conducted an RNA-seq analysis. The transcriptomic analysis revealed that the upregulation of E4 Type esterase gene NI-EST1 contributed to the elevated resistance of BPH19 in comparison to BPH80 and BPH15. These findings suggest that NI-EST1 plays a pivotal role in the increased resistance of N. lugens to carbamate insecticides fenobucarb.