Skip to main content

A common mechanism of detoxification for lambda-cyhalothrin and abamectin in Cydia pomonella

Abstract

Background

The primary method utilized by orchard owners to combat Cydia pomonella is the application of various chemical insecticides. However, this has resulted in the development of resistance. The resistance mechanisms to insecticides from different chemical classes are diverse but interconnected. Therefore, it is crucial to comprehend the commonalities in these mechanisms to effectively develop strategies for managing resistance.

Materials and methods

To determine whether target-site insensitivity to LCT and AM plays a role in resistance, the sequences of voltage-gated sodium channel (VGSC) and glutamate-gated chloride channel (GluCl) containing the mutation domains were detected. To validate whether similar mechanisms were involved in the detoxification process of lambda-cyhalothrin (LCT) and abamectin (AM) at sublethal doses (specifically LD10 and LD30), cytochrome P450 monooxygenases (P450), glutathione S-transferases (GST), and carboxylesterases (CarE) activities were evaluated after insecticides exposure; synergistic experiments were conducted using piperonyl butoxide (PBO), diethyl maleate (DEM), and triphenyl phosphate (TPP) as inhibitors of P450, GST, and CarE respectively. RNA-sequencing (RNA-Seq) was performed to compare the expression levels of detoxification-related genes between susceptible (SS) and resistant strains.

Results

The best known target-site mutations caused by LCT and AM, including L1014F in VGSC and V263I, A309V, I321T, and G326E in GluCl were not occurred. We observed that PBO had a strong synergistic effect on LCT and AM, while DEM on LCT. The activities of detoxification enzymes increased after insecticide exposures, indicating that the detoxification of LCT was primarily carried out by P450 and CarE enzymes, while P450 and GST enzymes played a major role in the detoxification of AM. A total of 72 P450 genes and 75 CarE genes were identified in the C. pomonella transcriptome, with 43 of these genes (including 11 P450, 3 GST, 10 CarE, 11 ABC transporters, and 8 UDP-glycosyl transferases) being over-expressed in response to both insecticides Interestingly, ABC transporters were predominantly induced by AM treatment, while GST showed higher induction levels with LCT treatment. Furthermore, LCT-resistant strains of C. pomonella exhibited higher levels of induction of detoxification-related genes compared to susceptible strains.

Conclusion

The up-regulation of these detoxification genes is a common metabolic mechanism employed by C. pomonella to counteract the effects of insecticides, although the extent of gene expression change varies depending on the specific insecticide.

Background

The codling moth, Cydia pomonella (Lepidoptera: Tortricidae), is a globally notorious pest with a wide distribution (Reyes et al. 2007; Wang et al. 2019). The larvae of C. pomonella cause fruit abscission by tunnelling into the centre of fruits to feed on the sarcocarp and seeds, leaving frass on the surface (Reyes et al. 2007; Voudouris et al. 2011). Although environmentally-sound management methods have been developed to suppress C. pomonella in certain regions (Witzgall et al. 2008; Knight et al. 2022), their effectiveness against high-density pest populations remains often limited (Calkins and Faust, 2003). Due to the easy availability and immediate visible efficacy of chemical insecticides, orchard owners predominantly rely on chemical applications as their primary strategy, leading to the development of C. pomonella resistance to insecticides (Reyes et al. 2007; Rodríguez et al. 2011; Bosch et al. 2018). To date, resistance evolution in field populations of C. pomonella towards pyrethroids has been documented in at least eight countries (Yang and Zhang, 2015; Ju et al. 2021). For instance, a field population collected from an orchard in the USA has developed resistance (~ tenfold) to lambda-cyhalothrin (Mota-Sanchez et al. 2008), while three field populations from Argentina have exhibited high resistance (> 30-fold) to lambda-cyhalothrin (Soleño et al. 2020). Avermectins derived from macrocyclic lactones (MLs) have been extensively utilized for controlling C. pomonella (Ju et al. 2021). However, the increasing usage inevitably exerts selection pressure for developing C. pomonella resistance (Reyes et al. 2007; Rodríguez et al. 2011; Bosch et al. 2018).

Insects have developed multiple mechanisms to survive in the presence of insecticides, including target-site insensitivity to insecticides or enhanced metabolic detoxification (Li et al. 2007; Yang et al. 2016; Yang et al. 2017; Paula et al. 2021; Zhou et al. 2023). The former is caused by point mutations of target genes that lead to the development of insecticide resistance. Generally, pyrethroids target voltage-gated sodium channel (VGSC) and avermectins target glutamate-gated chloride channel (GluCl) of insects (Mermans et al. 2017; Ju et al. 2021). Compared to target resistance, metabolic detoxification is achieved through the overexpression of enzymes or multi-drug transporter molecules (Yang et al. 2016; Yang et al. 2017; Li et al. 2022a). Previous studies have confirmed that resistance to insecticides is mainly conferred by three enzymatic complexes, including the cytochrome P450 monooxygenases (P450), glutathione S-transferases (GST), and carboxylesterases (CarE) (Wheelock et al. 2005; Yang 2016; Yang et al. 2017; Wang et al. 2019; Hu et al. 2020a, b). Our laboratory’s previous research found that resistance to LCT in C. pomonella field population was associated with enhanced P450 and GST activity and reduced CarE activity (Wei et al. 2020; Li et al. 2023). Additionally, UDP-glycosyl transferases (UGT) and ATP-binding cassette (ABC) transporters also contribute to the detoxification processes by conjugating bio-transformed or native toxicants and mediating the efflux of toxic substances (Bock 2016; Jin et al. 2019; Ju et al. 2022). However, most previous research has focused on the detoxification or resistance mechanism(s) of insects to a certain insecticide or a specific group of insecticides. For instance, GST contributes to the detoxification of pyrethroids (Hu et al. 2020a, b, 2023), and ABC is involved in avermectins metabolism (Jin et al. 2019; Ju et al. 2022). Nevertheless, insects may encounter various insecticides with different modes of action, resulting in resistance development against multiple mechanisms (Li et al. 2022a). Therefore, it is crucial to elucidate the commonness of resistance mechanisms towards different types of insecticides and facilitate strategies for managing resistance.

Synergists are a valuable tool for investigating metabolic resistance because they can specifically inhibit the detoxification pathways of insecticides (Gonzalez-Morales and Romero 2019). Piperonyl butoxide (PBO), a p450 inhibitor, is a commonly used synergist that has been employed for decades to study metabolic resistance in arthropods and enhance the effectiveness of insecticides, particularly when issues of resistance arise (Assis et al. 2018; Gonzalez-Morales and Romero 2019). Other synergists such as triphenyl phosphate (TPP), which inhibits carboxylesterases, and diethyl maleate (DEM), which inhibits GST, are also utilized to determine whether metabolic resistance plays a role in arthropod resistance to insecticides (Lilly et al. 2016; Gonzalez-Morales and Romero 2019). Previous investigations have shown that the synergistic efficacy of specific inhibitors varies among different insecticides. For example, DEM is more effective when co-formulated with pyrethroids while PBO exhibits broader compatibility with various insecticide classes like pyrethroids and macrocyclic lactones (Gonzalez-Morales and Romero 2019). TPP has demonstrated significant synergy with deltamethrin in some insects but not in others (Chigure et al. 2018; Gonzalez-Morales and Romero 2019). Furthermore, the synergism between inhibitors and insecticides also varies among different species of insects (Tian et al. 2018; Gonzalez-Morales and Romero 2019; Zhao et al. 2020; Li et al. 2022b). However, the detoxification mechanisms for insecticides from different chemical classes remain poorly understood, and employing synergistic studies may serve as a powerful approach to explore the primary types of mechanisms underlying insecticide resistance.

Traditional insecticides toxicity assessments have focused on acute lethal effects, but recent research emphasizes the importance of considering sublethal effects (Desneux et al. 2007). Thus, understanding both acute and sublethal effects is crucial for integrated pest management and pesticide registration procedures. In this study, we detected targeted site (VGSC and GluCl) mutations in an LCT resistant strain to determine if target-site insensitivity play a major role in C. pomonella resistance. We then measured the activity of P450, GST, and CarE in larvae response to insecticides exposure, and evaluated the synergistic effect of PBO, DEM, and TPP with LCT and AM in a susceptible strain. Subsequently, we performed a transcriptome analysis on third-instar larvae of C. pomonella treated with sublethal doses (LD10 and LD30) of LCT and AM to obtain comprehensive expression profiles of detoxifying enzyme genes. All P450 and CarE genes were systematically identified and characterized as part of this process. To investigate induction mechanisms for detoxification enzyme genes in response to different insecticides, we comprehensively profiled transcriptional abundance patterns across various treatments/tissues associated with different types of insecticide exposure. Furthermore, we compared expression levels between LCT-resistant strains and susceptible strains of C.pomonella for further insights into the roles played by detoxification genes in insecticide resistance. Our findings contribute towards unravelling the metabolic mechanism underlying broad-spectrum insecticide resistance at the molecular level while also aiding in developing scientific and sustainable pest management programs (IPM).

Materials and methods

Insect and chemicals

Two C. pomonella strains were used in this study. The first is a susceptible strain (SS), which has been maintained in the laboratory in a growth chamber (MLR-352H-PC, Panasonic) over 60 generations without exposure to any insecticides. The rearing conditions are 26 ± 1 ℃, 60 ± 5% relative humidity, and a photoperiod of 16:8 h (L:D). The rearing method was described by Wang et al (2019). The second strain, ZW_R, is a field population of C. pomonella collected from a pear orchard in Zhangwu County, Liaoning Province, northeast China (122.08°E, 42.28°N). This strain has developed a moderate level (16.97-fold) of resistance to LCT (Wei et al. 2020). The ZW_R strain was reared in the laboratory using the same conditions as the SS strain from July 2014 to October 2019. ZW_R strain was used for detection of target gene mutation, and expression analysis of detoxification genes in resistant population.The technical grade (> 98% purity) lambda-cyhalothrin and abamectin was purchased from Aladdin Reagent (Shanghai, China). The synergists used piperonyl butoxide (PBO, > 99% purity) as a P450 inhibitor, diethyl maleate (DEM, 98% purity) as a GST inhibitor, and triphenyl phosphate (TPP, 98% purity) as a CarE inhibitor purchased from Aladdin Reagent (Shanghai, China).

Detection of point mutations of target genes

The mRNA sequences of GluCl from Plutella xylostella (GenBank: GQ221939.1) was queried against the C. pomonella genome (http://v2.insect-genome.com/Organism/224) and transcriptome (SRR14923516) to identify CpGluCl gene (Additional file 1: Figure S1). Protein sequences of insects GluCl transmembrane domain 1 and 3 (TM1 and TM3) were aligned using ClustalW.

The segments of the legs of the adults were used for DNA extraction. DNA extraction was conducted following the manufacturer's instruction (Tiangen, Beijing, China). Gene-specific primers of CpGluCl gene which containing V263, A309, I321, and G326E domains (Additional file 1: Figure S1), were designed by Primer Premier 6. The primer pair of CpVGSC containing L1014 domain was obtained from previous (Soleño et al. 2020). Primers sequence were shown in Additional file 1: Table S1.

PCR was conducted with a total reaction volume of 20 µL. The reaction mixture consisted of 10 μL of TaKaRa Ex Taq (TaKaRa, Dalian, China), 1 µL of each primer (1 µM), 1 µL of DNA template, and 7 µL of ddH2O. The qPCR cycling conditions were as follows: an initial denaturation at 98 ℃ for 2 min, followed by 35 cycles of denaturation at 98 ℃ for 30 s, annealing at 60 ℃ for 30 s, and extend at 72 ℃ for 1 min, finally extend at 72 ℃ for 10 min to obtain PCR products. After sequencing, the nucleotide sequences of ZW_R and SS were compared to detect any target site mutations.

Bioassay

The bioassay was conducted using the topical application method with an Eppendorf pipettor (Hamburg, Germany) as described by Rodríguez et al. (2011). Third-instar larvae of the SS strain were used for this study. The lambda-cyhalothrin was dissolved in acetone and then diluted to eight concentrations (0, 4, 8, 16, 32, 64, 128, and 256 mg L−1). A volume of 2 μL of each concentration of insecticide solution was applied on the thoracic dorsum of third instar larvae. Controls were treated only with acetone following the same procedure. Fifteen larvae were treated at each concentration, and each treatment was repeated three times. The treated larvae were individually reared in a 2 × 5 cm glass tube and fed with artificial diet. Mortality was recorded at 12 h, 24 h, and 48 h post-exposure (hpe). A larva was considered dead if it did not respond to stimulation by an ink brush. The mortality data was underwent Probit analyses (Abbott 1925). The slopes and lethal dose (LD) values were estimated according to Finney (1971) using the Ld-p Line® software.

The bioassay of abamectin has been previously documented by Ju et al. (2022). For the exposure experiment, the LD10, LD30, and LD50 values were used, which were determined to be 4.657, 19.348, and 65.640 ng larvae−1 at 72 h, respectively.

Synergistic effects of detoxification enzyme inhibitors

To conduct synergism tests in SS strain, PBO, DEM, and TPP were dissolved in acetone to create a stock solution with a concentration of 300 mg/L, 300 mg/L, and 1000 mg/L, respectively. A 100 μL of this stock solution was then applied on the surface of a 0.5 × 0.5 × 0.5 cm artificial diet in each well of a 24-well plate. After allowing the solution to air-dry for 1 h in a fume hood, third instar larvae, starved for 12 h, were placed individually individually on the pretreated artificial diet. A control group was also included, consisting of larvae treated with an equal amount of acetone solution. One hour after the synergists treatments, each insecticide (LCT or AM) was applied onto the thoracic dorsum of each third-instar larva to conduct the bioassay as previously described. Mortality was recorded at 24, 48, and 72 h post-exposure, and the results at 48 h post-exposure for LCT or 72 h post-exposure for AM were used for the analysis of the synergism effect. The synergism ratio (SR) was calculated as the ratio of the lethal dose 50 (LD50) of the larvae treated with the test insecticide to the LD50 of the larvae treated with the test insecticide plus the synergists.

Sublethal insecticide exposure

In the exposure experiment, we utilized sublethal doses of LCT and AM, specifically LD10 and LD30, which are doses that would result in the mortality of 10% and 30% of the larvae, respectively. To administer the insecticide solution, we carefully applied two μL of each concentration on the thoracic dorsum of third-instar larvae. Acetone-treated larvae were used as controls and underwent the same procedure. Increase in weight of l3 larvae was measured after each sublethal insecticides treatments by averaging the weight of random 10 larvae as a group, at different time points: 24, 48, and 72 h after exposure. Additionally, surviving larvae were collected and frozen in liquid nitrogen for subsequent experiments at 12, 24, 48, and 72 hpe. Each treatment consisted of three replicates, with 10 larvae for each time interval.

For the analysis of tissue-specific expression, LD30 of LCT and AM was administered on third-instar larvae. At the 48-h mark, we dissected the head (HE), cuticle (CU), fat body (FB), midgut (MD), and Malpighian tubules (MT) from three groups of 30 larvae, and stored them at – 80 ℃.

Determination of detoxifying enzymes activity

Five l3 larvae were selected from each treatments with 3 repeats for the preparation of enzyme sources, the buffers for P450, GST, and CarE were 100 mM sodium phosphate (pH 7.8), 50 mM sodium phosphate (pH 7.2) and Tris–HCl/CaCl2 (25 mM/ 1 mM, pH7.0, containing 1 mM dithithreitol), respectively. BCA Protein Assay Kit (TaKaRa, Dalian, China) was used to determine the concentration of enzyme solution.

The activity of P450 was determined according to Tang et al. (2011), using 7-ethoxycoumarin as a substrate to construct a standard curve, and the specific activity of ECOD was calculated. CDNB was used as a substrate to determine GST activity according to the methods of Li et al. (2023). CarE activity was measured using a Tiangen Kit (Beijing, China) following the manufacturer’s instructions.

RNA isolation and RNA-Seq

Total RNA was extracted from entire body of larvae treated with sublethal LCT or AM as described in 2.5, three replicates, each pooled five larvae. The TaKaRa MiniBEST Universal RNA Extraction Kit (TaKaRa, Dalian, China) was employed, following the manufacturer’s instructions. The concentration and integrity of the RNA samples were assessed using a NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), respectively. RNA samples with an OD260/280 value ≥ 1.8 and integrity value ≥ 7.0 were considered suitable for further steps. These qualified RNA samples were then utilized for mRNA preparation and cDNA library construction. Subsequently, NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) was used to generate the sequencing libraries, following the manufacturer’s recommendations. The quality of the libraries was assessed using the Agilent Bioanalyzer 2100 system. Before analyzing the data, raw data of fastq format were firstly processed through in-house perl scripts. In this step, clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30 and GC content the clean data were calculated. All downstream analyses were based on clean data with high quality. The reference genome and gene model annotation files were directly downloaded from the genome website. The index of the reference genome was built using Hisat2 v2.0.5, and the paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5 (http://ccb.jhu.edu/software/hisat2/index.shtml) (Kim et al. 2015). Differential expression analysis between the two groups was conducted using the DESeq2 R package (1.20.0). DESeq2 utilizes statistical methods based on the negative binomial distribution to determine differential expression in digital gene expression data. The resulting P-values were adjusted using Benjamini and Hochberg's approach to control the false discovery rate. Differentially expressed genes (DEGs) were identified based on a Padj < 0.05 and |log2 (Fold Change)|> 1 cut-off for further analysis.

Identification of detoxifying enzyme genes from C. pomonella

The amino acid sequences of P450 and CarE from Bombyx mori (Xia et al. 2004; Yu et al. 2009) were obtained from the NCBI GenBank (http://www.ncbi.nlm.nih.gov/genbank/). These sequences were then queried against the C. pomonella genome and transcriptome database (SRR14923516) using the local BLASTP program. The E-value threshold of 10–5 was used to identify any potential P450 and CarE genes.

The location of the C. pomonella P450 and CarE genes on the chromosome was obtained from the InsectBase (http://www.insect-genome.com/). A phylogenetic tree was constructed using the maximum likelihood (ML) phylogenetic inference with IQ-TREE, integrated into the software PhyloSuite using an Edge-linked partition model with a bootstrap value of 5000 replicates (Minh et al. 2013; Nguyen et al. 2015).

cDNA synthesis and RT-qPCR analysis

The synthesis of cDNA was performed by utilizing 1 µg of RNA from each treatment, employing the PrimeScript RT reagent Kit with gDNA Eraser (Takara Dalian, China). The resulting cDNA was stored at a temperature of − 20 ℃ until it was ready for use. Gene-specific primers (Additional file 1: Table S1) were designed for each detoxification gene using Primer Premier 6. RT-qPCR was conducted on a Bio-Rad CFX96 (BioRad, U.S.A.) with a total reaction volume of 20 µL. The reaction mixture consisted of 10 μL of TB Green Premix Ex Taq 2 (TliRNaseH Plus; TaKaRa, Dalian, China), 1 µL of each primer (1 µM), 1 µL of cDNA template, and 7 µL of ddH2O. The RT-qPCR cycling conditions were as follows: an initial denaturation at 95 ℃ for 30 s, followed by 40 cycles of denaturation at 95 ℃ for 5 s, and annealing at 60 ℃ for 30 s. A melting curve analysis was performed starting at 63 ℃ for 5 s, with increments of 0.5 ℃ up to 95 ℃ to assess the specificity of the PCR products. The EF-1α (MN037793) and β-Actin (KC832921) of C. pomonella were used as reference genes (Wei et al. 2020). Each sample was subjected to three technical replicates and three biological replicates. The relative expression levels of each detoxification gene were determined using the 2CT method (Livak and Schmittgen 2001). Statistical significance was assessed using the one-way analysis of variance (ANOVA) with Duncan's test, using SPSS Statistics 20 software (IBM corporation, Armonk, NY, USA). All data are reported as the mean of triplicates ± standard deviation (SD). The significance level was set at P < 0.05.

Result

Detection of the target mutation in LCT-resistant population

Sequence alignment of GluCl subunits revealed that position V263 which located in TM1, and A309, I321, and G326 located in TM3 were highly conserved in all lepidoptera species (Additional file 1: Figure S2).

Sequencing results showed that in the CpGluCl gene, the codons of GTA, GCC, ATC and GGC encode valine (V) at 263, alanine (A) at 309, isoleucine (I) at 321 and glycine (G) at 323, respectively. Nucleotide mutations were not found at any of these sites (Additional file 1: Figure S3A). Likewise, leucine (L) at position 1014 of the VGSC gene is encoded by TTT, and there is no nucleotide mutation at this site (Additional file 1: Figure S3B).

Toxicity of LCT and AM against C. pomonella larvae

The linear regression of the dose-mortality relationship in SS larvae was fitted for LCT, and the Probit analyses are listed in Table 1. The LD10, LD30, and LD50 values were 4.724, 10.596, and 18.541 ng larvae−1 at 48 h, respectively (Table 1). The mortality in the controls was lower than 10%, suggesting that the bioassay results were valid. As for AM, the linear regression of the dose-mortality relationship has been documented by Ju et al. (2022).

Table 1 Dose–Response data for C. pomonella to lambda-cyhalothrin over Different Times

Relationship between detoxifying enzymes and insecticides metabolism in C. pomonella larvae

Detoxifying enzymes activity response to inecticides stress

The enzyme activity of P450 was up-regulated after two doses of two insecticides, and a time-dependent increase were observed (Fig. 1A). As for GST, enzyme activity increased after different LCT treatments and reached the highest level at 24 hpe in both two dose treatments; while GST activity was only increased significantly in LD30 of AM treatment at 72 hpe (Fig. 1B). In addition, CarE activity decreased first and then increased after both two insecticides treatments (Fig. 1C). It is worth noting that compared to the two dose treatments, the change in detoxification enzyme activities after LD10 treatment were higher than that of LD30 treatment.

Fig. 1
figure 1

The results of enzymatic assays. A P450 activities (B) GST activities (C) CarE activities. Blue histogram represents control; green histogram represent lambda-cyhalothrin treatment; pink histogram represent abamectin treatment. Letters on the error bars indicate significant differences between treatments by ANOVA analysis (P < 0.05)

Synergists increased the sensitivity of C. pomonella larvae to insecticides

We investigated the toxicity of LCT or AM in combination with the synergists PBO, DEM, and TPP, respectively. The results revealed that pre-treatment with the GST inhibitor DEM and the cytochrome P450 inhibitor PBO significantly increased the toxicity of LCT in the SS population, with synergistic ratio (SR) values of 3.793 and 3.139, respectively. However, the synergistic effect of TPP on LCT toxicity was less pronounced compared to other two synergists, with an SR of 1.455 (Table 2).

Table 2 Dose–Response data for C. pomonella to lambda-cyhalothrin and abamectin without and with piperonyl butoxide (PBO), diethyl maleate (DEM) and triphenyl phosphate (TPP)

For AM, the LD50 value was reduced from 65.640 ng larvae−1 to 15.258 ng larvae−1 for third-instar larvae of C. pomonella when AM combined with PBO, resulting in the highest SR of 4.302 among all treatments. Our result also indicated that the TPP synergist enhanced the toxicity of AM in larvae, with an SR value of 1.616. The combination of GST inhibitor DEM showed the lowest SR value among all tested synergists with AM.

Effects of LCT and AM on the body mass of C. pomonella

The body mass of the codling moth was measured after exposure to two sublethal doses of LCT and AM for different durations (24 h, 48 h, and 72 h). The weight gain of third instar larvae exhibited a significant decrease in LCT and AM exposure compared to the control (CK), except for the AM-LD10 treatment at 24 h. The impact of the LD10 dose of LCT on weight gain was greater than that of the LD30 dose, but only at 24 h did it show a significant difference (Fig. 2A).

Fig. 2
figure 2

Summary of the weight and number of differentially expressed genes in third-instar C. pomonella larvae were examined after exposure to LD10 and LD30 of LCT or AM at various time intervals. A Increase in weight of l3 in each treatment. Blue represents control; the hollow green histogram and the solid green histogram represent LD10 and LD30 of lambda-cyhalothrin treatment, respectively; hollow pink histogram and the solid pink histogram represent LD10 and LD30 of abamectin treatment, respectively. B The number of up-regulated and down-regulated genes in each treatment. Green and pink represent up-regulated genes and down-regulated genes compared to the control. LCT, lambda-cyhalothrin; AM, abamectin

Transcriptome response to LCT and AM treatment

To investigate the potential response mechanism of C. pomonella to sublethal doses of LCT and AM stress, RNA-seq was conducted on the codling moth treated with the two sublethal doses for different time periods (12 h, 24 h, and 48 h for LCT; 24 h, 48 h, and 72 h for AM). A total of 48 samples, including insecticide-treated and control samples, yielded an average of 64,775,370 raw sequencing reads, of which 63,640,539 clean reads were obtained after filtering low-quality reads (Additional file 1: Table S2). The clean reads ratio was consistently above 97.53%, with an average of 97.99%, meeting base call quality at Q20, indicating good quality of the clean reads (Additional file 1: Table S2). The GC% ranged from 41.59% to 48.43%, with an average of 41.54% for all samples.

The histogram provides an overview of the differentially expressed genes (DEGs) between the treatment groups and the control group for two insecticides (Fig. 2B). The results demonstrate a time-dependent increase in gene induction in C. pomonella larvae after exposure to LCT. Specifically, the number of up-regulated DEGs increases from 7900 at 24 h to 9319 at 48 h and 9811 at 72 h after exposure to LD10 of LCT. Similarly, the number of up-regulated DEGs increases from 7489 at 24 h to 11,059 at 48 h and 13,161 at 72 h after exposure to LD30 of LCT. Conversely, the number of down-regulated DEGs decreases over time. In the case of AM treatment, DEGs tend to be induced rather than inhibited with both sublethal doses at different time intervals, particularly at 48 h. At this time point, there are 12,855 and 12,284 up-regulated DEGs at LD10 and LD30 treatments, respectively, compared to only 9102 and 9516 down-regulated genes.

Furthermore, a total of 3168 DEGs are co-upregulated in LD10 LCT treatment and 1477 DEGs are co-upregulated in LD30 LCT at all time points (12 h, 24 h, and 48 h). Among these, 598 genes are co-upregulated in both sublethal doses of LCT. GO enrichment analysis reveals that the DEGs are primarily enriched in biological process (BP), cellular component (CC), and molecular function (MF), with a particular emphasis on biological process (Fig. 3A). Additionally, to further assess the effects of LCT on C. pomonella, the down-regulated DEGs at different time points (12, 24, and 48 h) after exposure are also counted. A total of 5916 and 2361 DEGs are down-regulated in LD10 and LD30 LCT treatment, respectively. The Venn diagram in Additional file 1: Figure S4A shows that there are 1131 genes that are co-downregulated in the LCT treatment.

Fig. 3
figure 3

The most enriched Gene ontology of co-upregulated expressed genes in the LCT-treated or AM-treated group. A In each section of the Venn diagram is the number of co-upregulated genes for LD10 and LD30 of LCT treatment at different time intervals, compared to the control. Gene ontology (GO) enriched DEGs of C. pomonella after exposure to two sublethal doses of LCT at all different time intervals. The Y-axis represents sub-Go terms of biological process (BP), cellular component (CC), and molecular function (MF). The significance of pathways is shown as q-value (colour bar), rich factor (X-axis), and circle size representing the number of genes regulated. B In each section of the Venn diagram is the number of co-upregulated genes for LD10 and LD30 of AM treatment at different time intervals, compared to the control. Gene ontology (GO) enriched DEGs of C. pomonella after exposure to two sublethal doses of AM at all different time intervals. The Y-axis represents sub-Go terms of biological process (BP), cellular component (CC), and molecular function (MF). The significance of pathways is shown in q-value (colour bar), rich factor (X-axis), and circle size, showing the number of genes regulated. LCT, lambda-cyhalothrin; AM, abamectin

After treatment with AM, the number of up-regulated differentially expressed genes (DEGs) was higher than the number of down-regulated DEGs, which is in contrast to LCT. LD10 and LD30 AM exposure resulted in the induction of 4228 DEGs and 3372 DEGs, respectively, compared to the control. Venn diagrams revealed that 1734 genes were co-upregulated. The majority of DEGs were associated with biological processes and molecular function, while only a few DEGs were related to cellular components (Fig. 3B). Additionally, the number of down-regulated DEGs in response to LD10 AM treatment was higher than that of LD30, with 2965 and 2460 DEGs, respectively. Common down-regulated DEGs (1236) were observed after treatment with both LD10 and LD30 AM, and gene ontology (GO) analyses demonstrated significant enrichment in molecular function (Additional file 1: Figure S4B).

The fold change response to the two insecticides ranged from 12.85 to 14.89, and the top 10 up-regulated DEGs in each treatment, as measured by fold change, are listed in Additional file 1: Table S3. The majority of the top 10 up-regulated DEGs in most treatments were related to insect cuticle protein, trypsin, and GMC oxidoreductase. Additionally, several genes involved in detoxification processes were strongly induced. For example, the CarE gene, CPOM18541, exhibited a 6.41-fold and 8.37-fold increase after 24 and 48 h of LD10 LCT exposure, respectively. CPOM13402 showed an 11.20-fold and 11.06-fold change after exposure to both LD10 and LD30 AM for 24 h, respectively. Furthermore, the transcript abundance of the UDP-glucosyl transferase gene (CPOM21547) was up-regulated by over fivefold in LCT-LD10-12 h, LCT-LD10-24 h, and LCT-LD30-24 h treatments compared to the control.

Identification and classification of detoxifying enzyme families

P450 gene family

Following an extensive search and analysis, we have successfully identified a total of 72 P450 genes from C. pomonella. These genes have been categorized into four distinct Clans, specifically referred to as Clan 2, Clan 3, Clan 4, and mitochondrial Clan. A thorough phylogenetic analysis has revealed that P450 genes belonging to the same clan form tight clusters. Among these clusters, Clan 3 exhibited the highest number of genes with 31, followed by Clan 4 with 26 genes and mitochondrial Clan with 9 genes. Lastly, Clan 2 comprised the smallest cluster, featuring only 6 genes (Additional files 1: Figure S5, Table S4). Furthermore, we conducted an analysis to determine the distribution of P450 genes on the genomes’ scaffolds. The results showed that these 72 genes were spread across 26 different chromosomes (Additional file 1: Figure S6).

CarE gene family

We identified 75 CarE genes in C. pomonella through local BLAST searches. Subsequent NCBI blast and phylogenetic analysis allowed us to classify these CarEs into three classes, each containing multiple clades. The intracellular catalytic class (clade: α-esterase) comprised the largest number of genes with 53, while the secreted catalytic class (clade: juvenile hormone esterase, integument esterase, β-esterase, and uncharacterized) and neurodevelopmental class (clade: acetylcholinesterase, uncharacterized, gliotactin, neuroligin, and neurotactin) contained 11 genes each (Additional files 1: Figure S7, Table S5). The 75 CarE genes were found on a total of 20 different chromosomes, with chromosomes 4, 8, and 16 hosting the largest clusters of CarE genes (Additional file 1: Figure S8).

Expression analysis of detoxification genes

Expression analysis of detoxification genes in response to LCT and AM exposure

To investigate the potential contribution of detoxification enzyme genes in the metabolism of two insecticides, we conducted a screening process on 598 and 1734 co-expressed genes after exposing the subjects to sublethal doses (LD10 and LD30) of LCT and AM at different time points, and 9 and 13 detoxifying enzyme genes were found (Fig. 3). In addition, 3 highly induced genes were screened out in Additional file 1: Table S3, so a total of 25 detoxifying genes were selected for further analysis. The expression patterns of these genes in insecticide-treated larvae at various time points are illustrated in Fig. 4A. Our findings indicate that the 9 induced genes responding to LCT did not exhibit significant up-regulation (< 1.5-fold) after AM treatment, and in fact, some of them were down-regulated. However, the 13 genes that showed consistent up-regulation in response to AM treatment were also induced by LD30 of LCT treatment at 48 h, with the exception of CPOM14343 (CarE) which decreased by 0.27-fold. Moreover, both LCT and AM treatments led to increased transcriptional abundance of the 3 highly induced genes (Fig. 4A).

Fig. 4
figure 4

The expression profiles of selected detoxification enzyme genes were co-upregulated and expressed in the LCT-treated or AM-treated group. A Transcriptomic analysis (RNA-seq) of detoxification enzyme genes caused by LCT and AM, log2(fold change) were shown on a heatmap with colours ranging from green (low expression) to pink (high expression). Treatments are shown above, gene family information and gene id were given on the left. Genes with an underline mean they are selected from the LCT-treated group; genes without an underline mean they are selected from the AM-treated group; genes with a grey background mean that they were selected from the top 10 up-regulated unigenes in LCT or AM treatments. B In each section of the Venn diagram is the number of co-upregulated genes for LD10 or LD30 of LCT or AM treatment at 48 hpe, compared to the control. C Transcriptomic analysis (RNA-seq) of detoxification enzyme genes screened by Venn diagram ( B) caused by LCT and AM, log2(fold change) were showed on heatmap with colours ranging from white (low expression) to pink (high expression). Treatments are shown above, gene family information and gene id were given on the left

As a result of the induction of most up-regulated DEGs by both insecticides at 48 h (Figs. 2B, 4A), we developed Veen to investigate the common detoxification mechanisms at the molecular level (Fig. 4B). Out of the 4793 co-upregulated genes, a total of 21 genes are associated with detoxification metabolism and have an induction level of more than onefold. Three of these genes (CPOM18541 (CarE), CPOM02212 (CarE), and CPOM11913 (UGT)) overlap with the previously screened genes and are not repeated in Fig. 4C to avoid duplication. Figure 4C displays the transcriptome data of the remaining 18 co-upregulated detoxification genes. The majority of these genes exhibited increased transcriptional abundance in response to the stress caused by the two insecticides. Furthermore, genes that were highly induced after LCT treatment also showed a strong increase after AM treatment, including the P450 gene (CPOM02891, CPOM15879), CarE gene (CPOM128050), and UGT gene (NOVEL.3146), with expression levels being nearly fourfold higher than the control or even more. Additionally, the GST genes (CpGSTe7 and CpGSTe8) were more strongly induced by LCT compared to AM (Fig. 4C).

Expression analysis of detoxification genes in tissues

The expression patterns of the genes that were screened in various tissues (head, cuticle, fat body, midgut, and Malpighian tubule) of larvae were analyzed using RT-qPCR after being exposed to sublethal doses (LD10 and LD30) of LCT and AM for 48 h (Fig. 5). The results indicated that the transcript abundance of these genes showed similar profiles across the tested tissues, with lower expression levels observed in the cuticle and higher expression levels in the fat body. However, four P450 genes (CPOM09449, CPOM19073, CPOM14496, and CPOM08165) exhibited higher expression in the head compared to other tissues after exposure to both insecticides. Most P450 genes showed higher expression levels in the fat body, midgut, and Malpighian tubules following exposure to LCT and AM. Two GST genes, CpGSTe7 and CpGSTe8, were predominantly induced in the fat body (12.56- and 12.67-fold, respectively) when third-instar larvae were exposed to LD10 of LCT and AM. Another GST gene, CpGSTe4, showed high induction in the midgut compared to other tissues after exposure to LCT. The expression levels of several detoxification enzyme genes, including CPOM06460 (UGT), CPOM00266 (UGT), CPOM06856 (UGT), and CPOM03338 (ABC), were strongly induced in the fat body, midgut, and Malpighian tubule when third-instar larvae were exposed to both LD10 and LD30 of LCT and AM. Notably, regardless of the insecticide treatment, LD10 treatment generally had a greater impact on the transcript abundance level between tissues compared to LD30 treatment, showing a more pronounced up-regulation level (Fig. 5).

Fig. 5
figure 5

Tissue-specific expression profiles of detoxification enzyme genes in the third-instar larvae of C. pomonella after exposure to LD10 and LD30 of LCT or AM for 48 h measured by RT-qPCR, log2(fold change) with colours ranging from green (low expression) to pink (high expression). Treatments are given above: Head, HE; Cuticle, CU; Fat body, FB; Midgut, MD; Malpighian tubules, MT. Gene family information and gene id were given on the left. Genes with underlined means they are selected from the LCT-treated group; genes without underlined means they are selected from the AM-treated group; genes with a grey background mean that they were selected from the top 10 up-regulated unigenes in LCT or AM treatments, genes in red font means they were co-upregulated genes for LD10 or LD30 of LCT and AM treatment at 48 hpe, compared to the control

Expression analysis of detoxification genes in the resistant population

To determine if these detoxification genes played a role in insecticide resistance, their transcript abundance was analyzed in different C. pomonella strains with varying levels of resistance to LCT. The analysis revealed that, compared to the SS strain, the expression levels of CPOM11145 and CPOM09449 genes were lower in the resistance strain. However, the transcript abundance of other tested P450 genes was significantly higher in the ZW_R strain, with most of them showing over a twofold induction. Similarly, the transcriptional levels of CpGSTe4 and CpGSTe7 genes were 6.94-fold and 5.84-fold increased in the ZW_R strain compared to SS, respectively. The CarE family genes also exhibited a similar pattern, with 6 out of 10 genes showing over a 4.92-fold increase in expression in the ZW_R stain compared to SS. Additionally, the expression level of 8 ABC genes in the ZW_R population was significantly higher than in SS, with the ABCC subfamily gene, CPOM07642, showing an 18.64-fold higher transcript abundance in the ZW_R population. The UGT genes also showed a similar expression pattern, with the increase in transcript abundance being correlated with the resistance levels (Fig. 6).

Fig. 6
figure 6

The expression level of detoxification enzyme genes in different C. pomonella strains measured by RT-qPCR. Genes with an underline mean they are selected from the LCT-treated group; genes without an underline mean they are selected from the AM-treated group; genes with a grey background mean that they were selected from the top 10 up-regulated unigenes in LCT or AM treatments, genes in red font means they were co-upregulated genes for LD10 or LD30 of LCT and AM treatment at 48 hpe, compared to the control. Green bars represent susceptible strain, and pink bars represent ZW_R strain, the data were normalized to the expression of EF-1α and β-Actin, and the change in the expression level was calculated using the 2−ΔΔCt method

Discussion

Insect resistance to insecticides poses a complex challenge to agricultural productivity. An understanding of how insects develop resistance to these chemicals, as well as their potential fitness costs, is crucial in order to delay this process, manage resistance in insects, and prolong the effectiveness of insecticides (Amezian et al. 2021; Saeed et al. 2021; Gui et al. 2022; Ju et al. 2022; Mocchetti et al. 2023). By uncovering the common resistance mechanisms employed by insects against insecticides with diverse modes of action, it becomes possible to establish a strategy for managing pest populations that have developed resistance to multiple insecticides simultaneously.

Reduced insecticide efficacy in C. pomonella has been related to non-specific mechanisms including enhanced enzymatic metabolization and/or modification of the specific molecular target of one group of insecticides (Soleño et al. 2020; Ju et al. 2021). The increased activity of P450 and GST and the decreased activity of CarE has been suggested to explain resistance to LCT in ZW_R strain, in line with prevoius researches in C. pomonella (Li et al. 2023). In this study, we did not detect any mutations in the target genes of two insecticides in ZW_R strain, the disappearance or/and lack of mutations may be a consequence due to the absence of ABM selection pressure or a high VGSC mutation-introduced fitness cost in C. pomonella (Ni et al 2023). Thus, we propose that enzymatic-mediated metabolic resistance is major mechanism of low level resistance in C. pomonella strain.

Synergistic effects and changes in detoxification enzyme activity often contribute to specify the major metabolic mechanisms in various insects (Gonzalez-Morales and Romero 2019; Sun et al. 2017; Shen et al. 2020). In this particular study, we not only found that the activity of P450 enzyme in SS larvae increased after both LCT and AM stress, but also found PBO, a potent inhibitors of P450, significantly increased the toxicity of LCT (SR = 3.793) and AM (SR = 4.302). This finding is consistent with previous research that has shown the involvement of P450 in resistance development against pyrethroids and macrocyclic lactones, as evidenced by synergism studies with PBO (Tian et al. 2018; Zhao et al. 2020). Although TPP was not highly synergistic toxicity to insecticides, the increased CarE activity was also associated with sensitivity of SS larvae to LCT and AM, which was inconsistent with CarE decrease in ZW_R (Tian et al. 2018; Zhao et al. 2020; Li et al. 2023). This may be related to the differences in sensitivity between strains, and the mechanism of CarE involvement in insecticide resistance is also varies between increasing or decreasing (Ju et al. 2021). GST contributes to the detoxification of LCT were clearly proved in both resisitance and sensitive strains before, our study also confirmed that LCT exposure inducted the GST activity of larvae, while DEM also had the highest synergistic effects (SR = 4.302) on LCT (Wang et al. 2019; Hu et al. 2020a, b). Additionally, synergism between DEM and ivermectin had also been reported in a few reasearches, while we only found a weak synergistic effect when DEM combined with AM and a few induction on GST activity after AM exposure (Wang et al. 2016; Shakya et al. 2022).

The transcriptome results from our study demonstrate a time-dependent pattern of gene induction in C. pomonella larvae following exposure to LCT and AM. Specifically, there was an increase in the number of up-regulated genes from 12 to 48 h for LCT, and from 24 to 48 h for AM. These findings align with previous research (Wei et al. 2019), indicating that the upregulation of detoxifying genes is a common mechanism in insects to regulate intracellular insecticide concentration. Additionally, we observed that the number of down-regulated genes was higher in the LD30 treatment compared to the LD10 treatment. Moreover, exposure to LD30 dose of insecticide had a more negative impact on weight gain than LD10 dose. This aligns with the findings of Khan et al. (2021), who suggested that the higher number of down-regulated differentially expressed genes (DEGs) may contribute to adverse effects on body weight, which was confirmed in our study. These results indicate that exposure to insecticides induces significant molecular and biological changes in C. pomonella larvae, and increased gene expression may serve as a mechanism for insects to respond to external insecticides.

The GO annotations of commonly up-regulated DEGs, in response to LCT exposure, indicate a preference for enrichment in biological processes, particularly metabolic processes. This suggests that these genes may be involved in the metabolism of LCT (Wei et al. 2019; Amezian et al. 2021). On the other hand, the commonly up-regulated DEGs after AM exposure show enrichment not only in metabolic processes but also in cell communication and signal signaling processes. This finding demonstrates that there are both similarities and differences in the detoxification response between LCT and AM in C. pomonella. The reason for this difference lies in the main mechanism of action of AM, which involves inhibiting the transmission of electrical impulses in the muscles and nerves of invertebrates by amplifying the effects of glutamate on the invertebrates-specific gated chloride channel (Bloomquist 2003). Consequently, after AM exposure, larvae are more likely to regulate cell communication and signal signaling processes by altering the abundance of gene expression.

To reveal the detoxification mechanism of C. pomonella in response too LCT and AM, we focused on the analysis of three major detoxification enzyme genes. Among them, P450 and CarE have not been identified systematically in C. pomonella. In this study, a total of 72 P450s were identified from C. pomonella (Additional file 1: Table S4). The number of genes is comparable to other related species such as Bombyx mori (78) (Xia et al. 2004), Spodoptera litura (68) (Hu et al. 2019), and Grapholita Molesta (77) (Guo et al. 2017), but less than Tribolium Castaneum (143) (Zhu et al. 2013) and Anopheles gambiae (105) (Pavlidi et al. 2013) (Additional file 1: Table S4). This is mainly because B. mori, S. litura, and G. Molesta all belong to Lepidoptera, of which G. molesta and C. pomonella both belong to Noctidae and have closely relationship shared similar ecological niches. Previous studies have indicated that Clan 3 and Clan 4 are the most diverse, while the smaller Clan 2 and mitochondrial Clan are considered ortholog genes across species (Feyereisen et al. 2012). In our study, we observed a greater number of gene families in the CYP3 and CYP4 Clans compared to the other two Clans, and gene duplication and divergence were found to be the main driver of CYPome evolution. Furthermore, these genes in the two Clans have been shown to play a positive role in the adaptive xenobiotic metabolism and detoxification of insecticides (Nauen et al. 2022). As for CarEs, we identified 75 genes divided into three classes. Among them, the catalytic class contains the highest number of genes (53), which is similar to other Lepidoptera species such as B. mori (55) (Yu et al. 2009). Many of the closely catalytic class genes are part of lineage-specific family expansions and are arrayed in clusters on chromosomes, including chr4, chr8, and chr16. The expansion of catalytic class genes has been linked to the degradation of plant volatiles or other xenobiotics, including insecticides (Wheelock et al. 2005; Yu et al. 2009; Nauen et al. 2022). Our findings suggest that the subfamilies associated with detoxification metabolism consist of a large number of genes that cluster together, indicating their potential involvement in insecticide detoxification in C. pomonella. Additionally, we identified fewer P450 and CarE genes compared to the genome annotation results of Tong et al (2022). This discrepancy may be attributed to the consideration of genome assemblies as parameters only, without filtering, leading to a bias in the gene family size. Therefore, we believe that our results are more comprehensive and accurate.

Previous studies have confirmed that the induction of the detoxification genes plays a crucial role in adaption to various insecticides stresses, including but not limited to LCT and AM (Wheelock et al. 2005; Amezian et al. 2021; Ju et al. 2022; Messenger et al. 2021; Nauen et al. 2022). Our transcriptome results demonstrate that gene expression in each detoxification enzyme family is upregulated after exposure to both LCT and AM, suggesting that the tolerance of C. pomonella to insecticides is not regulated by a single gene, but rather by multiple gene regulation. Detoxification genes are involved in responding to multiple insecticides, not just a single one. Specifically, ABC genes account for a higher proportion of commonly up-regulated detoxification genes in all doses and time points of AM treatments. On the other hand, three common upregulated GST genes were observed to have higher upregulation after LCT treatments compared to AM treatments. It has been reported that increased expression of ABC genes is a mechanism observed in insects, including C. pomonella and Helicoverpa armigera, to regulate the concentration of insecticides inside their cells, particularly abamectin (Jin et al. 2019; Ju et al. 2022). Additionally, previous studies have shown that the expression of GST genes in C. pomonella is significantly upregulated under LCT and is involved in the phase II detoxification of pyrethroids by conjugating the thiol group of glutathione (GSH) to molecules possess an electrophilic centre (Yang and Zhang, 2015; Wang et al. 2019; Hu et al. 2020a,b; Ju et al. 2021). These findings indicate that the upregulation of detoxification genes is a common mechanism for detoxification, but the number of upregulated genes varies depending on the insecticide.

In order to better understand the role of detoxification genes in C. pomonella’s ability to metabolize insecticides, we conducted an investigation into the tissue-specific expression profiles, with a focus on metabolic tissues. Our findings revealed that the majority of these genes exhibit significant overexpression in the fat body, midgut, and Malpighian tubules when exposed to sublethal doses of LCT stress. The fat body serves as the primary organ for detoxification, while the midgut possesses the largest detoxification and metabolic tissues. Additionally, the Malpighian tubules play a role in the excretion of xenobiotics. This information, supported by previous studies (Lei et al. 2014; Mao et al. 2019; Hu et al. 2020a, b), strongly suggests that the overexpression of detoxification enzyme genes in these specific tissues serves as a broad-spectrum defense mechanism against LCT and AM in C. pomonella.

It is necessary to verify whether screened detoxification genes play a crucial role in regulating the expression level to adapt and overcome insecticides in the insecticide-resistant strain of C. pomonella. Previous studies have reported the overexpression of multiple P450s, GSTs, and CarEs in multi-insecticide-resistant populations of Anopheles arabiensis, which have been shown to metabolize various combinations of type I and type II pyrethroids (Messenger et al. 2021). Similarly, investigations have demonstrated the expression of numerous GST and UGT genes in Paederus fuscipes after exposure to a sublethal concentration of emamectin benzoate (Khan et al. 2021). In the genus Spodoptera insects, the superfamilies of detoxification genes, including P450, GST, UGT, and ABC, have been found to play a crucial role in adapting to insecticides and have been verified in different classes of insecticide-resistant strains. Among these detoxification genes, cytochrome P450 has been identified as the most upregulated, as reported in the literature (Amezian et al. 2021). Previous research has also documented that up-regulated genes in multi-insecticide resistance are involved in resistance (Amezian et al. 2021; Tong et al 2022), which aligns with our findings. Specifically, our study observed the overexpression of 8 P450 genes, 2 GST genes, 6 CarE genes, 8 ABC genes, and 6 UGT genes in an LCT resistance strain (ZW_R, RF:16.97). These results provide molecular evidence supporting the role of detoxification gene expression in insecticide resistance in C. pomonella.

Collectively, despite the distinct mechanisms of action exhibited by the two insecticides, they both exhibit a shared detoxification metabolic pattern, albeit with certain variations depending on the specific insecticide. Consequently, it is imperative to consider the potential cross-resistance between these two insecticides when formulating integrated pest management strategies. Furthermore, it can be postulated that C. pomonella populations that have developed resistance to one insecticide are more prone to developing resistance upon exposure to other classes of insecticides; a factor that should not be underestimated.

Conclusion

The findings of our study demonstrate the significant synergistic effect of PBO in combating both LCT and AM, as well as the effectiveness of DEM against LCT and TPP against LCT and AM in C. pomonella. Besides, increased detoxification enzymes activity was observed in C. pomonella response to insecticide exposures, this provides valuable and directional insight into the crucial role of detoxification enzymes in the breakdown of LCT and AM. Additionally, our study successfully identified the P450s and CarEs from the genome database and transcriptomic of C. pomonella for the first time, systematically shedding light on the detoxification genes involved. Our results suggest that these detoxification gene superfamilies play a significant role in increasing transcript abundance and facilitating the broad-spectrum detoxification of various insecticides. It is important to consider cross-resistance in the management and control of C. pomonella based on our findings.

Availability of data and materials

The authors are willing to share all data and additional information on materials used in this study upon a written request to the corresponding author.

References

  • Amezian D, Nauen R, Le Goff G. Comparative analysis of the detoxification gene inventory of four major Spodoptera pest species in response to xenobiotics. Insect Biochem Mol Biol. 2021;138:103646.

    Article  CAS  PubMed  Google Scholar 

  • Assis CP, Gondim MG, Siqueira HA. Synergism to acaricides in resistant Neoseiulus californicus (Acari: phytoseiidae), a predator of Tetranychus urticae (Acari: tetranychidae). Crop Prot. 2018;106:139–45.

    Article  CAS  Google Scholar 

  • Bloomquist JR. Chloride channels as tools for developing selective insecticides. Arch Insect Biochem Physiol. 2003;54(4):145–56.

    Article  CAS  PubMed  Google Scholar 

  • Bock KW. The UDP-glycosyltransferase (UGT) superfamily expressed in humans, insects and plants: animal-plant arms-race and co-evolution. Biochem Pharmacol. 2016;99:11–7.

    Article  CAS  PubMed  Google Scholar 

  • Bosch D, Rodríguez MA, Depalo L, Avilla J. Determination of the baseline susceptibility of European populations of Cydia pomonella (Lepidoptera: Tortricidae) to chlorantraniliprole and the role of cytochrome P450 monooxygenases. J Econ Entomol. 2018;111(2):844–52.

    Article  CAS  PubMed  Google Scholar 

  • Calkins CO, Faust RJ. Overview of areawide programs and the program for suppression of codling moth in the western USA directed by the United States Department of Agriculture-Agricultural Research Service. Pest Manag Sci. 2003;59(6–7):601–604.

  • Chigure GM, Sharma AK, Kumar S, Fular A, Sagar SV, Nagar G, Upadhaya D, Saravanan BC, Kumar R, Ghosh S. Role of metabolic enzymes in conferring resistance to synthetic pyrethroids, organophosphates, and phenylpyrazole compounds in Rhipicephalus microplus. Int J Acarol. 2018;44:28–34.

    Article  Google Scholar 

  • Desneux N, Decourtye A, Delpuech JM. The sublethal effects of pesticides on beneficial arthropods. Ann Rev Entomol. 2007;52(1):81-106.ju.

    Article  CAS  Google Scholar 

  • Feyereisen R. Insect CYP genes and P450 enzymes insect molecular biology and biochemistry. Cambridge: Academic Press; 2012. p. 236–316.

    Book  Google Scholar 

  • Gonzalez-Morales MA, Romero A. Effect of synergists on deltamethrin resistance in the common bed bug (Hemiptera: Cimicidae). J Econ Entomol. 2019;112(2):786–91.

    Article  CAS  PubMed  Google Scholar 

  • Gui F, Lan T, Zhao Y, Guo W, et al. Genomic and transcriptomic analysis unveils population evolution and development of pesticide resistance in fall armyworm Spodoptera frugiperda. Protein Cell. 2022;13(7):513–31.

    Article  CAS  PubMed  Google Scholar 

  • Guo Y, Chai Y, Zhang L, Zhao Z, Gao LL, Ma R. Transcriptome analysis and identification of major detoxification gene families and insecticide targets in Grapholita Molesta (Busck) (Lepidoptera: Tortricidae). J Insect Sci. 2017;17(2):43.

    Article  PubMed  PubMed Central  Google Scholar 

  • Hu B, Zhang SH, Ren MM, Tian XR, Wei Q, Mburu DK, Su JY. The expression of Spodoptera exigua P450 and UGT genes: tissue specificity and response to insecticides. Insect Sci. 2019;26(2):199–216.

    Article  CAS  PubMed  Google Scholar 

  • Hu C, Wang W, Ju D, Chen GM, Tan XL, Mota-Sanchez D, Yang XQ. Functional characterization of a novel lambda-cyhalothrin metabolising glutathione S-transferase, CpGSTe3, from the codling moth Cydia pomonella. Pest Manag Sci. 2020a;76(3):1039–47.

    Article  CAS  PubMed  Google Scholar 

  • Hu C, Wei ZH, Li PR, Harwood JD, Li XY, Yang XQ. Identification and functional characterization of a sigma glutathione s-transferase CpGSTs2 involved in lambda-cyhalothrin resistance in the codling moth Cydia pomonella. J Agric Food Chem. 2020b;68(45):12585–94.

    Article  CAS  PubMed  Google Scholar 

  • Hu C, Liu YX, Zhang SP, Wang YQ, Gao P, Li YT, Yang XQ. Transcription factor AhR regulates glutathione S-transferases conferring resistance to lambda-Cyhalothrin in Cydia pomonella. J Agric Food Chem. 2023;71(13):5230–9.

    Article  CAS  PubMed  Google Scholar 

  • Jin M, Liao C, Chakrabarty S, Zheng W, Wu K, Xiao Y. Transcriptional response of ATP-binding cassette (ABC) transporters to insecticides in the cotton bollworm, Helicoverpa armigera. Pestic Biochem Physiol. 2019;154:46–59.

    Article  CAS  PubMed  Google Scholar 

  • Ju D, Mota-Sanchez D, Fuentes-Contreras E, Ya-Lin Z, Xiao-Qi W, Xue-Qing Y. Insecticide resistance in the Cydia pomonella (L): Global status, mechanisms, and research direction. Pestic Biochem Physiol. 2021;9:104925.

    Article  Google Scholar 

  • Ju D, Dewer Y, Zhang S, Hu C, Li P, Yang X. Genome-wide identification, characterization, and expression profiling of ATP-binding cassette (ABC) transporter genes potentially associated with abamectin detoxification in Cydia pomonella. Ecotoxicol Environ Saf. 2022;230:113152.

    Article  CAS  PubMed  Google Scholar 

  • Khan MM, Khan AH, Ali MW, Hafeez M, Ali S, Du C, Fan Z, Sattar M, Hua H. Emamectin benzoate induced enzymatic and transcriptional alternation in detoxification mechanism of predatory beetle Paederus fuscipes (Coleoptera: Staphylinidae) at the sublethal concentration. Ecotoxicology. 2021;30(6):1227–41.

    Article  CAS  PubMed  Google Scholar 

  • Kim D, Landmead B, Salzberg SL. Hisat: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357-U121.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  • Knight AL, Preti M, Basoalto E, Mujica MV, Favaro R, Angeli S. Combining female removal with mating disruption for management of Cydia pomonella in apple. Entomologia Generalis. 2022;42:309–21.

    Article  Google Scholar 

  • Li XC, Schuler MA, Berenbaum MR. Molecular mechanisms of metabolic resistance to synthetic and natural xenobiotics. Annu Rev Entomol. 2007;52:231–53.

    Article  PubMed  Google Scholar 

  • Li W, Wang X, Jiang P, Yang M, Li Z, Huang C, He Y. A full-length transcriptome and gene expression analysis of three detoxification gene families in a predatory stink bug picromerus lewisi. Front Physiol. 2022a;13:1016582.

    Article  PubMed  PubMed Central  Google Scholar 

  • Li R, Zhu B, Liang P, Gao XW. Identification of carboxylesterase genes contributing to multi-insecticide resistance in Plutella xylostella (L.). Entomol Generalis. 2022b;42:967–76.

    Article  Google Scholar 

  • Li PR, Shi Y, Ju D, Liu YX, Wang W, He YS, Zhang YY, Yang XQ. Metabolic functional redundancy of the CYP9A subfamily members leads to P450-mediated lambda-cyhalothrin resistance in Cydia pomonella. Pest Manag Sci. 2023;79(4):1452–66.

    Article  CAS  PubMed  Google Scholar 

  • Lei Y, Zhu X, Xie W, Wu Q, Wang S, Guo Z, Xu B, Li X, Zhou X, Zhang Y. Midgut transcriptome response to a Cry toxin in the diamondback moth, Plutella xylostella (Lepidoptera: Plutellidae). Gene. 2014;533(1):180–187.

  • Lilly DG, Dang K, Webb CE, Doggett SL. Evidence for metabolic pyrethroid resistance in the common bed bug (Hemiptera: Cimicidae). J Econ Entomol. 2016;109(3):1364–8.

    Article  CAS  PubMed  Google Scholar 

  • Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25(4):402–8.

    Article  CAS  PubMed  Google Scholar 

  • Mao T, Li F, Fang Y, Wang H, Chen J, Li M, Li Z, Qu J, Li J, Hu J, Chen X, Ni M, Li B. Effects of chlorantraniliprole exposure on detoxification enzyme activities and detoxification-related gene expression in the fat body of the silkworm, Bombyx mori. Ecotoxicol Environ Saf. 2019;176:58–63.

    Article  CAS  PubMed  Google Scholar 

  • Mermans C, Dermauw W, Geibel S, Van Leeuwen T. A G326E substitution in the glutamate-gated chloride channel 3 (GluCl3) of the two-spotted spider mite Tetranychus urticae abolishes the agonistic activity of macrocyclic lactones. Pest Manag Sci. 2017;73(12):2413–2418.

  • Messenger LA, Impoinvil LM, Derilus D, Yewhalaw D, Irish S, Lenhart A. A whole transcriptomic approach provides novel insights into the molecular basis of organophosphate and pyrethroid resistance in Anopheles arabiensis from Ethiopia. Insect Biochem Mol Biol. 2021;139:103655.

    Article  CAS  PubMed  Google Scholar 

  • Minh BQ, Nguyen M, Haeseler AV. Ultrafast approximation for phylogenetic bootstrap. Mol Biol Evol. 2013;30(5):1188–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  • Mocchetti A, Dermauw W, Van Leeuwen T. Incidence and molecular mechanisms of insecticide resistance in Frankliniella occidentalis, Thrips tabaci and other economically important thrips species. Entomol Generalis. 2023. https://doi.org/10.1127/entomologia/2023/1889.

    Article  Google Scholar 

  • Mota-Sanchez D, Wise JC, Poppen RV, Gut LJ. Hollingworth R M. Resistance of codling moth, Cydia pomonella (L.) (Lepidoptera: Tortricidae), larvae in Michigan to insecticides with different modes of action and the impact on field residual activity. Pest Management Science, 2008;64(9):881–890.

  • Nauen R, Bass C, Feyereisen R, Vontas J. The role of cytochrome P450s in insect toxicology and resistance. Annu Rev Entomol. 2022;67:105–24.

    Article  CAS  PubMed  Google Scholar 

  • Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32(1):268–74.

    Article  CAS  PubMed  Google Scholar 

  • Ni R, Wang Y, Zhong Q, Li M, Zhang D, Zhang Y, Qiu X. Absence of known knockdown resistance mutations but fixation of CYP337B3 was detected in field populations of Helicoverpa armigera across China. Pestic Biochem Physiol. 2023;195:105542.

    Article  CAS  PubMed  Google Scholar 

  • Paula DP, Lozano RE, Menger JP, Andow DA, Koch RL. Identification of point mutations related to pyrethroid resistance in voltage-gated sodium channel genes in Aphis glycines. Entomol Generalis. 2021;41:243–55.

    Article  Google Scholar 

  • Pavlidi N, Dermauw W, Rombauts S, Chrisargiris A, Leeuwen TV, Vontas J. Analysis of the olive fruit fly Bactrocera oleae transcriptome and phylogenetic classification of the major detoxification gene families. PLoS ONE. 2013;8:e66533.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  • Reyes M, Franck P, Charmillot PJ, Ioriatti C, Olivares J, Pasqualini E, Sauphanor B. Diversity of insecticide resistance mechanisms and spectrum in European populations of the codling moth Cydia Pomonella. Pest Manag Sci. 2007;63(9):890–902.

    Article  CAS  PubMed  Google Scholar 

  • Rodríguez MA, Bosch D, Avilla J. Resistance of Spanish codling moth (Cydia pomonella) populations to insecticides and activity of detoxifying enzymatic systems. Entomol Exp Appl. 2011;138(3):184–92.

    Article  Google Scholar 

  • Saeed R, Abbas N, Hafez AM. Biological fitness costs in emamectin benzoate-resistant strains of Dysdercus koenigii. Entomol Generalis. 2021;41:267–78.

    Article  Google Scholar 

  • Shakya M, Nandi A, Fular A, Kumar S, Bisht N, Sharma AK, Singh K, Kumar R, Kumar S, Juliet S, Ghosh S. Synergistic property of piperonyl butoxide, diethyl maleate, triphenyl phosphate and verapamil hydrochloride with deltamethrin and ivermectin against Rhipicephalus microplus ticks. Ticks Tick-Borne Dis. 2022;13(6):102006.

    Article  PubMed  Google Scholar 

  • Shen J, Li Z, Li D, Wang R, Zhang S, You H, Li J. Biochemical mechanisms, cross-resistance and stability of resistance to metaflumizone in Plutella xylostella. Insects. 2020;11(311):11050311.

    Google Scholar 

  • Soleño J, Parra-Morales LB, Cichón L, Garrido SA, Montagna CM. Occurrence of pyrethroid resistance mutation in Cydia pomonella (Lepidoptera: Tortricidae) throughout Argentina. Bull Entomol Res. 2020;110(2):201–6.

  • Sun H, Pu J, Chen F, Wang J, Han Z. Multiple ATP-binding cassette transporters are involved in insecticide resistance in the small brown planthopper, Laodelphax striatellus. Insect Mol Biol. 2017;26(3):343–355

  • Tang BZ, Sun JY, Zhou XG, Gao XW, Liang P. The stability and biochemical basis of fufenozide resistance in a laboratory-selected strain of Plutella xylostella. Pestic Biochem Physiol. 2011;101:80–5.

    Article  CAS  Google Scholar 

  • Tian F, Mo X, Rizvi SAH, Li C, Zeng X. Detection and biochemical characterization of insecticide resistance in field populations of Asian citrus psyllid in Guangdong of China. Sci Rep. 2018;8(1):1–11.

    Article  Google Scholar 

  • Tong D, Zhang L, Wu N, et al. The oriental armyworm genome yields insights into the long-distance migration of noctuid moths. Cell Rep. 2022;41(12):111843.

    Article  CAS  PubMed  Google Scholar 

  • Voudouris CC, Sauphanor B, Franck P, Reyes M, Mamuris Z, Tsitsipis JA, Vontas J, Margaritopoulos JT. Insecticide resistance status of the codling moth Cydia pomonella (Lepidoptera: Tortricidae) from Greece. Pestic Biochem Physiol. 2011;100(3):229–38.

    Article  CAS  Google Scholar 

  • Wang X, Wang R, Yang Y, Wu S, O’Reilly AO, Wu Y. A point mutation in the glutamate-gated chloride channel of Plutella xylostella is associated with resistance to abamectin. Insect Mol Biol. 2016;25(2):116–25.

    Article  CAS  PubMed  Google Scholar 

  • Wang W, Hu C, Li XR, Wang XQ, Yang XQ. CpGSTd3 is a lambda-Cyhalothrin metabolizing glutathione S-transferase from Cydia pomonella (L.). J Agric Food Chem. 2019;67(4):1165–72.

    Article  CAS  PubMed  Google Scholar 

  • Wei N, Zhong Y, Lin L, Xie M, Zhang G, Su W, Li C, Chen H. Transcriptome analysis and identification of insecticide tolerance-related genes after exposure to insecticide in Sitobion avenae. Genes. 2019;10(12):951.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  • Wei ZH, Liu M, Hu C, Yang XQ. Overexpression of glutathione S-transferase genes in field lambda-Cyhalothrin-resistant population of Cydia pomonella: reference gene selection and expression analysis. J Agric Food Chem. 2020;68(21):5825–34.

    Article  CAS  PubMed  Google Scholar 

  • Wheelock CE, Shan G, Ottea J. Overview of carboxylesterases and their role in the metabolism of insecticides. J Pestic Sci. 2005;30(2):75–83.

    Article  CAS  Google Scholar 

  • Witzgall P, Stelinski L, Gut L, Thomson D. Codling moth management and chemical ecology. Annu Rev Entomol. 2008;53:503–22.

    Article  CAS  PubMed  Google Scholar 

  • Xia Q, Zhou Z, Lu C, Cheng D, Yang HM. A draft sequence for the genome of the domesticated silkworm (Bombyx mori). Science. 2004;306:1937–40.

    Article  PubMed  Google Scholar 

  • Yang XQ. Gene expression analysis and enzyme assay reveal a potential role of the carboxylesterase gene CpCE-1 from Cydia pomonella in detoxification of insecticides. Pestic Biochem Physiol. 2016;129:56–62.

    Article  CAS  PubMed  Google Scholar 

  • Yang XQ, Zhang YL. Characterization of glutathione S-transferases from Sus scrofa, Cydia pomonella and Triticum aestivum: Their responses to cantharidin. Enzym Microb Technol. 2015;69:1–9.

  • Yang XQ, Wang W, Tan XL, Wang XQ, Dong H. Comparative analysis of recombinant cytochrome P450 CYP9A61 from Cydia pomonella expressed in Escherichia coli and Pichia pastoris. J Agric Food Chem. 2017;65(11):2337–44.

    Article  CAS  PubMed  Google Scholar 

  • Yu QY, Lu C, Li WL, Xiang ZH, Zhang Z. Annotation and expression of carboxylesterases in the silkworm. Bombyx Mori BMC Genom. 2009;10(1):1–14.

    Google Scholar 

  • Zanger UM, Schwab M. Cytochrome P450 enzymes in drug metabolism: regulation of gene expression, enzyme activities, and impact of genetic variation. Pharmacol Ther. 2013;138:103–41.

    Article  CAS  PubMed  Google Scholar 

  • Zhao YX, Huang JM, Ni H, Guo D, Yang FX, Wang X, Wu S, Gao CF. Susceptibility of fall armyworm, Spodoptera frugiperda (JE Smith), to eight insecticides in China, with special reference to lambda-cyhalothrin. Pestic Biochem Physiol. 2020;168:104623.

    Article  CAS  PubMed  Google Scholar 

  • Zhou YX, Han Q, Feng K, Wang JY, Zhou HF, Wen M, Duan HX, Wang YL, Ren BZ. New binding sites of nicotinic acetylcholine receptors from Myzus persicae. Entomol Generalis. 2023. https://doi.org/10.1127/entomologia/2023/1662.

    Article  Google Scholar 

  • Zhu F, Moural TW, Shah K, Palli SR. Integrated analysis of cytochrome P450 gene superfamily in the red flour beetle. Tribol Castaneum BMC Genom. 2013;14:174.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This research was supported by the National Key R&D Program of China (2021YFD1400200), and the National Natural Science Foundation of China (32272588).

Author information

Authors and Affiliations

Authors

Corresponding authors

Correspondence to Yuting Li or Xueqing Yang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no conflict of interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

GluCl mRNA (A) and amino acid sequence (B) of C. pomonella, the specific primer positions were underlined in red, the PCR product was in red letter, positions of V263I, A309V, I321T, and G326E mutations were highlight with yellow, red, blue, and green backgrounds. Figure S2. Alignment of transmembrane domain 2 and 3 (TM2 and TM3) of C. pomonella and other species, including Plutella xylostella (PxGluCl-V263 and PxGluCl-263I, GenBank: ACT09139.1,), Bombyx mori (BmGluCl, AGW27406.1), Drosophila melanogaster (DmGluCl, AAG40735.1), Musca domestica (MdGluCl, BAD16657.1) and Tetranychus urticae (TuGluCl, NP_001310061.1). Figure S3. Nucleotide sequence of GluCl and VGSC of C. pomonella in SS and ZW_R. Figure S4. The most enriched Gene ontology (GO) of co-downregulated expressed genes in LCT-treated or AM-treated group. A In each section of the Venn diagram is the number of co-downregulated genes for LD10 and LD30 of LCT treatment at different time intervals, compared to the control. Gene GO enriched DEGs of C. pomonella after exposure to two low different doses of LCT at all different time intervals. The Y-axis represents sub GO terms of biological process (BP), cellular component (CC), and molecular function (MF). The significance of pathways is shown in q-value (color bar), rich factor (X-axis), and circle size, showing the number of genes regulated. B In each section of the Venn diagram is the number of co-downregulated genes for LD10 and LD30 of AM treatment at different time intervals, compared to the control. GO enriched DEGs of C. pomonella after exposure to two different doses of AM at all different time intervals. The Y-axis represents sub Go terms of biological process (BP), cellular component (CC), and molecular function (MF). The significance of pathways is shown in q-value (color bar), rich factor (X-axis), and circle size, showing the number of genes regulated. LCT lambda-cyhalothri; AM abamectin. Figure S5. Phylogenetic tree of P450 genes in Cydia pomonella. Amino acid sequences from Bombyx mori (Bm), Spodoptera litura (Sl), Spodoptera exigua (Se), Spodoptera frugiperda (Sf) Manduca sexta (Ms), Plutella xylostella (Px), Zygaena filipendulae (Zf) and Helicoverpa armigera (Ha) were aligned using ClustalW (http://www.mbio.ncsu.edu/bioedit/bioedit.html). Phylogenetic analysis was constructed using the maximum likelihood with IQ-TREE. Four P450 clans are clustered together and are indicated by different colors. Figure S6. Protein location of the 72 P450s in C. pomonella. Genes names and chromosome numbers are showed on the right and left of the bar, respectively. Figure S7. Phylogenetic tree of CarE genes in Cydia pomonella. Amino acid sequences from Bombyx mori (Bm) and Chilo suppressalis (Csup) were aligned using ClustalW (http://www.mbio.ncsu.edu/bioedit/bioedit.html). Phylogenetic analysis was constructed using the maximum likelihood with IQ-TREE. Three CarE clans are clustered together and are indicated by different colors. Figure S8. Protein location of the 75 CarEs in C. pomonella. Genes names and chromosome numbers are showed on the right and left of the bar, respectively. Table S1 Primer sequence used for PCR or RT-qPCR. Table S2 Summary statistics of the transcriptome of C. pomonella. Table S3 Top 10 up-regulated unigenes. Table S4 Comparison of the number of gene in P450 subfamilies from various insect species. Table S5 Comparison of the number of gene in CarE subfamilies from various insect species.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ju, D., Hu, C., Li, P. et al. A common mechanism of detoxification for lambda-cyhalothrin and abamectin in Cydia pomonella. CABI Agric Biosci 4, 52 (2023). https://doi.org/10.1186/s43170-023-00192-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s43170-023-00192-0

Keywords