Genome Mining and Expression Analysis of Carboxylesterase and Glutathione S-Transferase Genes Involved in Insecticide Resistance in Eggplant Shoot and Fruit Borer, Leucinodes orbonalis (Lepidoptera: Crambidae)

The shoot and fruit borer, Leucinodes orbonalis (Lepidoptera: Crambidae) is the major cause of low productivity in eggplant and insecticides being the mainstay of management of L. orbonalis. However, field control failures are widespread due to the evolution of insecticide resistance. Taking advantage of the whole genome sequence information, the present study investigated the level of insecticide resistance and the expression pattern of individual carboxylesterase (CE) and glutathione S-transferases (GSTs) genes in various field collected populations of L. orbonalis. Dose-mortality bioassays revealed a very high level of resistance development against fenvalerate (48.2–160-fold), phosalone (94-534.6-fold), emamectin benzoate (7.2–55-fold), thiodicarb (9.64–22.7-fold), flubendiamide (187.4–303.0-fold), and chlorantraniliprole (1.6–8.6-fold) in field populations as compared to laboratory-reared susceptible iso-female colony (Lo-S). Over-production of detoxification enzymes viz., CE and GST were evident upon enzyme assays. Mining of the draft genome of L. orbonalis yielded large number of genes potentially belonging to the CE and GST gene families with known history of insecticide resistance in other insects. Subsequent RT-qPCR studies on relative contribution of individual genes revealed over-expression of numerous GSTs and few CEs in field populations, indicating their possible involvement of metabolic enzymes in insecticide resistance. The genomic information will facilitate the development of novel resistance management strategies against this pest.


INTRODUCTION
Eggplant or aubergine (Solanum melongena), also called as brinjal, is one of the most popular vegetable crops in southeast Asia especially India, China, and Bangladesh. The most important limiting factor in brinjal cultivation is the damage caused by shoot and fruit borer, Leucinodes orbonalis Guenee (Lepidoptera: Crambidae). It was first described from India and now it is distributed all over Asia, Africa, and in few parts of Europe (Mally et al., 2015). In India and Bangladesh, it causes severe yield losses up to 93% despite best management practices (Kodandaram et al., 2017;Prodhan et al., 2018). Surveys in India and Bangladesh indicated that farmers spray chemical insecticides up to 84 times during a 6-7 month cropping season. A very high level of pesticide load at a concentration 40-450 times higher than the maximum residue limit (MRL) in marketed brinjal fruits poses major health concerns (Srinivasan, 2008;Latif et al., 2010;Kariyanna et al., 2020b). It is also considered as a pest of quarantine significance to the number of countries outside its native range. L. orbonalis has high reproductive potential with overlapping generations. Studies reporting high level of insecticide resistance development in L. orbonalis might be attributed to the long-term indiscriminate use of various insecticides on eggplant (Kodandaram et al., 2017;Shirale et al., 2017).
Resistant insects overcome the toxic effect of insecticide molecules by adopting one or more mechanisms ranging from cuticular thickening, nerve impenetration, altered production of metabolic enzymes, target site insensitivity, enhanced excretion through ABC transporters as well as by gut symbionts . The enhanced detoxification of insecticides by metabolic enzymes viz., carboxylesterases (CE) and glutathione S-transferases (GSTs) are commonly reported in resistant populations of many insect species and mites (Ranson et al., 2002;Enayati et al., 2005;Oakeshott et al., 2005;Perry et al., 2011;Das and Dutta, 2014;Akıner and Ekşi, 2015;Benoit et al., 2015). The resistant individuals metabolize the insecticides faster due to higher catalytic rate or enhanced production of the enzyme as a consequence of increased transcription or gene duplication (Panini et al., 2016;Kariyanna et al., 2020a). The GSTs are phase II metabolic enzymes and are distributed in most of the animals (Hayes and Pulford, 1995). Among several subclasses of GST, the enzymes which are frequently associated with the insecticide resistance are Delta and Epsilon classes (Friedman, 2011;Lumjuan et al., 2011). Similarly, enhanced metabolism of organophosphates, carbamates, and pyrethroids are frequently associated with CEs through gene amplification, upregulation, mutations by coding sequence, or a combination of all these mechanisms (Zhang et al., 2012;Cui et al., 2015).
Considering the economic and social impact of L. orbonalis, genetically modified eggplant expressing insecticidal cry1Ac gene from Bacillus thuringiensis has been developed (Shelton et al., 2018), but, a moratorium on its commercialization was imposed by Indian Government due to public concern. Though Bt eggplant is commercialized in Bangladesh way back in 2014, yet the area under cultivation is <2500 Ha (Shelton et al., 2018) due to various reasons. Baring Bt brinjal, the insecticides remain the sole method of controlling L. orbonalis in all the eggplant growing countries. With the availability of genome sequence, the present study investigated the level of insecticide and the expression pattern of CE and GST genes from field collected insecticide resistant L. orbonalis populations, to pinpoint the key metabolic genes involved in insecticide degradation.

Insect Collection and Maintenance
The field populations of L. orbonalis larvae were collected during 2017-2018 from intensive eggplant growing regions of India viz., Raichur (16.2120  All the field populations were reared under laboratory conditions at 27 ± 2 • C, 60-70% relative humidity (RH), and a photoperiod of 14:10 h (L:D) on a natural diet and the F 1 individuals were used for bioassay and biochemical studies. The insecticide susceptible iso-female line (National Accession number: NBAIR-IS-CRA-01A), designated as Lo-S (65th generation) was originally derived from L. orbonalis collected near Bengaluru (12.9716 • N, 77.5946 • E) and maintained at insect genomic resources laboratory at ICAR-NBAIR.

Insecticide Resistance Bioassays
Insecticides fenvalarate (20% EC), emamectin benzoate (5% EC), phosalone (35%), thiodicarb (75%), flubendiamide (20%), and chlorantraniliprole (18.5%) were selected for conducting dosemortality bioassays based on their usage history on brinjal by farmers. Filter paper residue assay (Cheng et al., 2010) with essential modifications was used for insecticide resistance bioassays on early second instar larvae of L. orbonalis. Five to seven appropriate concentrations of each of the selected insecticides were prepared based on the dose bracketing technique. The filter paper discs (4.5 cm diameter) were dipped in the appropriate dilutions and dried vertically under shade for 1 h. Then, the discs were placed individually in plastic containers (5 cm diameter) and 10 larvae were released in each container. Moistened filter paper discs without insecticide were used in control. After 24 h, the larvae were transferred to untreated natural diet. All the assays were replicated and repeated at least thrice on alternate days. The mortality of the larvae was assessed after 48 h. The pooled larval mortality data were subjected to probit analysis using the software POLO (LeOra, 1987) and the lethal concentration to kill 50% of the test larvae (LC 50 ) was calculated for each population. Resistance ratios (RRs) were calculated with the following formula: LC 50 of field population/LC 50 of insecticide susceptible Lo-S population.

Preparation of Midgut Homogenate
The starved late second instar larvae were used for the preparation of midgut homogenate. Midguts were dissected and homogenized with homogenization buffer (0.1 M sodium phosphate buffer pH 7.8 containing 1 mM each of DTT, EDTA, PTU, and PMSF). The content was centrifuged at 12,000 rpm for 20 min and the clear supernatant was used as a enzyme source for estimating the titers of carboxylesterase and GST. The total protein content of the preparations was assessed by Coomassie brilliant blue G-250 dye-binding method using bovine serum albumin (BSA) as the standard (Bradford, 1976).

Assays on Carboxylesterase and Glutathione S-Transferase
The activity of carboxylesterase was determined using the method described by van Asperen (1962) with essential modifications and α-naphthyl acetate as a substrate. The assay mixture in a 96 well microplate (iMark, Biorad microplate reader) consisted of an appropriate amount of midgut homogenate, 0.1 M sodium phosphate buffer pH 7.8 and 800 µl of 3 mM α-naphthyl acetate containing 0.3 mM eserine. The mixture was incubated at 30 • C for 30 min. Finally, 200 µl of 0.1% tetrazotized o-dianisidine (Fast blue B) in 3.5% sodium dodecyl sulfate was added and incubated for 20 min at room temperature in dark. The α-naphthol formation was measured at 590 nm and the enzyme activity was computed from α-naphthol standard curve. Glutathione S-transferase activity was determined using 1-chloro-2,4-dinitrobenzene (CDNB) as a substrate (Kao et al.,  1989). The assay mixture composed of 50 mM each of CDNB and L-glutathione reduced and 10 µl gut homogenate. The change in absorbance was measured at 340 nm for 5 min and the enzyme activity in terms of nmoles of CDNB conjugated per minute per mg protein using the molar extinction coefficient of 5.3 mM −1 (optimized for the path length: 0.55 cm).

Esterase Isozyme Studies
The native polyacrylamide gel electrophoresis (PAGE) with 10% resolving gel was performed as per Davis (1964) to visualize the esterase isozymes from the 2nd instar midgut larval homogenates of the L. orbonalis populations tested. Midgut homogenates with 20 µg protein were loaded onto the native PAGE and run at a constant voltage of 70 for 1.5 h. Gel was briefly stained with freshly prepared 0.05% (w/v) α-naphthyl acetate and 0.1% (w/v) fast blue B in 50 mM phosphate buffer pH 7.8.

Genome Mining and Phylogenetic Analysis
The draft assembled genome (NCBI Bioproject ID: PRJNA377400) of L. orbonalis was used as query for extracting the gene sequences of glutathione S-transferases (GSTs) and carboxylesterases (CEs). The hmm model for carboxylesterase and GST (CE-PF00135, PF02230, and GST-C PF00043, GST-N PF02798, GST-PF13417, PF14497, PF17171, PF17172) were obtained from PFAM database 1 . A total of 94 and 33 putative CE and GST sequences were mined using hmm search with HMMER3 software package 2 . Among them, the genes with known history of involvement in insecticide resistance in other insects were selected and further confirmed with NCBI BLASTX. The GST and CE protein sequences of Drosophila melanogaster was downloaded from its genome database (Marygold et al., 2013), GenBank 3 and Uniprot 4 . We combined these putative GST protein sequences with D. melanogaster's annotated protein sequences for each of the six subclasses of GST (Epsilon, Sigma, Omega, Delta, Zeta, and Theta). The sequences were multiple aligned using CLUSTALW and a maximum likelihood tree with 1000 bootstrap value was constructed using MEGA7 5 and TreeDyn online server 6 was used to visualize the tree. Similarly, 94 CE protein sequences under seven subclasses (α-esterases, gliotactin, glutactin, neurolignin, juvenile hormone esterases, acetylcholinesterases, and cricklet) were combined, multiple aligned using MAFFT online server using conserved domain feature and phylogenetic tree of CE genes from L. orbonalis using D. melanogaster as a reference were constructed with the help of CLC workbench (QIAGEN Bioinformatics 7 ) using maximum likelihood approach.

Quantitative Real-Time PCR
The late 2nd instar larvae were used for RNA extraction (ISOLATE II RNA mini kit) by following the manufacturer's guidelines (Bioline). Purity and concentration of total RNA were measured in a spectrophotometer (NanoDrop Lite, Thermo Fisher Scientific) and denatured agarose gel (Masek et al., 2005). RNA samples with an A260/A280 ratio ranging from 1.8 to 2.0 and A260/A230 ratio >2.0 were used for cDNA preparation. First strand of complementary DNA was synthesized from 4 µg of total RNA using Revert AID first strand cDNA synthesis kit (Thermo Fisher Scientific TM , Lithuania) following the suppliers guidelines and stored at −80 • C till further use. qPCR primers for 25 GST and 16 CE gene sequences were designed using Primer 3.0 software. The detailed information of the primers used in the current study is listed in Tables 1, 2.
To remove the primer dimer and confirm the efficiency, OligoEvaluator TM sequence analysis tool (accessed on July 2018 8 ) was used. Primer specificity analysis was observed by single peak in the melting curve and single band from the agarose gel of the all candidate genes of CE and GST ranged from 90.36 to 104.58% (Tables 1, 2). Quantitative real-time PCR (RT-qPCR) amplifications were performed in 20 µL reaction consisted of 10 µL 2×SYBR R Premix EX TaqTM II (Tli RNaseH Plus, TAKARA R , Japan), 1 µl cDNA and gene specific primer pair. The parameters used for RT-qPCR are: One cycle of 95 • C for 3 min; 35 cycles of 95 • C for 30 s, 53 • C for 45 s, and 72 • C for 1 min; a final cycle of 72 • C for 10 min using Roche 480II machine. All the samples including control (no template) and internal control were performed in triplicates. The PCR products were electrophoresed on 2.0% agarose gel in a 1.0×TAE buffer. Relative expression levels for the CE and GST genes were calculated by 2 − CT method. The 28SR3 (28S ribosomal protein S3 mitochondrial)  gene was used as an internal control to normalize the expression of target genes.

Insecticide Resistance Monitoring
The LC 50 values against fenvalerate, phosalone, thiodicarb, emamectin benzoate, flubendiamide and chlorantraniliprole indicated that there is a large shift in susceptibility of fieldpopulations of L. orbonalis as compared to the laboratory reared susceptible iso-female colony (Lo-S). The field collected populations exhibited 3.6-160-fold resistance against fenvalerate, 7.3-534.6-fold resistance against phosalone, 7.0-55.0-fold resistance against emamectin benzoate, 2.0-22.7fold resistance against thiodicarb, 29.5-303.0-fold resistance against flubendiamide, and 1.6-8.6-fold resistance against chlorantraniliprole. However, the susceptibility level of all the field populations and the Lo-S did not vary significantly against chlorantraniliprole based on the overlapping fiducial limit values. High level of resistance development against phosalone and flubendiamide was observed in L. orbonalis collected from Varanasi whereas other field populations, the RRs were nonsignificant as compared to Lo-S. L. orbonalis populations from Bhubaneswar showed a significantly very high levels of resistance against fenvalerate, phosalone and flubendiamide as compared to Lo-S. The population of L. orbonalis collected from high hills near Pune recorded less RR and the susceptibility status was almost on par with Lo-S population ( Table 3).

Activities of Detoxification Enzymes
Quantitative differences in the titer of GST and CE were determined using α-naphthyl acetate and 1-chloro-2,4dinitrobenzene (CDNB) as substrates, respectively. Significant level of elevated GST activities (4.1-8.9-fold) was observed in field collected L. orbonalis as compared to susceptible Lo-S population. However, the overproduction of CE has not pronounced much (1.5-3.2-fold), but there was a significant level of overproduction in two field populations based on Mann-Whiteney U-test. The results indicated that GST activity was altered more profoundly than the CE ( Table 4).
The qualitative differences in carboxylesterase isozymes were visualized under native PAGE after incubation of the gels in the α-naphthyl acetate as a substrate. Among the four major esterase activity bands (E 1 to E 4 ), the high molecular weight E 1 and E 2 bands were absent in Lo-S and Pune populations and faint in case of Varanasi population. However, the E 1 band was very prominent and intense in other field populations indicating its over-production nature (Figure 1).

Differential Expression of CE Genes
Ninety-four CE genes were identified from the genome and transcriptome sequences of L. orbonalis. They were classified under seven subfamilies viz., α-esterases, gliotactin, glutactin, neurolignin, juvenile hormone esterases, acetylcholinesterases, and cricklet like orthologs (Figure 2). The length of the gene sequences ranged from 225 to 6498 bp. Expression profiling of 16 genes (homologous to resistant genes in other insects with ≥90% query coverage) was examined from the mRNA samples derived from late second instar larvae of five field-collected L. orbonalis populations by RT-qPCR. The CE genes contig11486-0.8 and contig4653-0.41 were over-expressed more than 10-fold as compared to susceptible Lo-S population (Figures 3, 4).

Differential Expression of Glutathione S-Transferase Genes
Mining of genome and transcriptome data of L. orbonalis yielded 33 unigenes codings for GSTs. Lengths of the identified GST unigene sequences ranged from 370 to 5,048 bp. The phylogenetic analysis revealed that the predicted GST genes are represented under all six classes of GST viz., Epsilon (10 genes), Sigma (6 genes), Omega (5 genes), Delta (8 genes), Zeta (3 genes), and Theta (1 gene) (Figure 5). RT-qPCR analysis of 25 GST genes (homologous to other insects resistant genes) revealed that many of them were over-expressed in larvae collected from the field. The GST genes named contig2323-0.22, contig14202-0.24, contig3755-0.39, contig3537-0.3, contig199-0.91, contig3023-0.16, contig3291-0.7, contig4841-1.1, and contig2323-0.23 were FIGURE 4 | Expression profiles (fold changes over Lo-S population) of carboxylesterase genes across field collected L. orbonalis populations (P, Pune; D, Dharmapuri; R, Raichur; B, Bhubaneswar; V, Varanasi). The fold changes are indicated in different shades indicate significant difference as per Mann-Whiteney U-test p-value < 0.05. FIGURE 5 | Phylogenetic tree of 34 predicted L. orbonalis glutathione S-transferase proteins with D. melanogaster genes representative of the six characterized GST classes. Putative GST proteins in L. orbonalis were identified using hmmsearch against PFAM GST families (GST-C PF00043, GST-N PF02798), the Maximum Likelihood tree with 1000 bootstraps was constructed using MEGA7 and TreeDyn online server was used to visualize the tree. Number on each node is bootstrap values in percent. Scale = 1 amino acid substitution per site.
very highly expressed (>50-fold) in one or more resistant field-collected L. orbonalis populations (Figures 6, 7) over the susceptible Lo-S population.

DISCUSSION
Damage caused by L. orbonalis is the major limiting factor in realizing maximum productivity of eggplant in many Asian countries including India, China, and Bangladesh. The high reproductive potential and shorter life cycle of the pest coupled with the evolution of multiple insecticide resistance pose a major challenge in profitable cultivation of eggplant (Prodhan et al., 2018). The insecticide resistance levels detected in the field collected populations ranged from low (>10-fold) to moderate (10-100-fold) and high (>100-fold) against fenvalerate (48.2-160-fold), phosalone (94-534.58-fold), emamectin benzoate (7-55-fold), thiodicarb (9.64-22.67-fold), and flubendiamide (29.51-363.91-fold). The resistance levels reflect the occurrence of differential selection pressure associated with long term history of insecticide use against L. orbonalis. However, the response to chlorantraniliprole was similar among the resistant field collected populations with RRs below 10fold. The efficacy of chlorantraniliprole has been reported earlier (Kodandaram et al., 2013;Munje et al., 2015). The chlorantraniliprole and flubendiamide are anthranilic and phthalic diamides, respectively, are the newer insecticides with little or no cross-resistance among them and also with other classes of insecticides. They are the activators of the ryanodine receptor.
Insect employs various mechanisms to nullify the toxic effect of xenobiotic compounds. Insect metabolic enzymes play a major  Frontiers in Physiology | www.frontiersin.org role in detoxification of plant allelochemicals and pesticidal compounds in addition to other physiological roles. Insect pests have evolved large reservoirs of various detoxification enzymes. Glutathione S-transferases mediated detoxification can be direct or by the metabolism of secondary products generated from other detoxification enzymes. The insecticides are metabolized rapidly into readily excretable water-soluble metabolites by reductive dehydrochlorination or by conjugation with reduced glutathione (Pavlidi et al., 2018). The carboxylesterase is a phase-1 detoxification enzyme that mainly hydrolyses organophosphates, carbamates, and synthetic pyrethroids.
The quantitative mechanisms underlying metabolic resistance is characterized by over-production due to gene amplification or transcriptional upregulation Feyereisen, 2015;Pavlidi et al., 2018). A qualitative mechanism occurs as a result of changes in the enzymatic characteristics. The present investigation on the enhanced metabolism of insecticides in resistant populations of L. orbonalis resistance is very prominent as inferred by elevated GST (1.2-2.6-fold) and carboxylesterases (6.3-13.1-fold) titers. The involvement of CEs and GSTs in insecticide resistance is commonly reported in many other insect species Pavlidi et al., 2018).
With the advent of genome and transcriptome information, multiple GSTs and carboxylesterases encoding genes have been characterized from insects and mites such as Plutella xylostella, Culex spp., Tribolium castaneum, Bombyx mori, Leptinotarsa decemlineata, Tetranychus cinnabarinus, Musca domestica and in many other insect species (Yang and Liu, 2011;Wang et al., 2014;Cui et al., 2015;Meisel and Scott, 2018;Pavlidi et al., 2018). L. orbonalis genome has large expansion of CE and GST genes which is comparable to many other insects. However, the GSTs have prominent role in insecticide resistance as compared to CEs. mRNA level of at least nine GSTs were over-expressed >50-fold where as none of the CEs showed very high levels of expression in the field collected resistant populations. The over-expression of few CE genes in field populations, especially from Varanasi, Raichur, and Bhubaneswar, also aligned with additional high molecular weight esterase band (E1) probably contributed to overall resistance. A similar high molecular weight esterase band associated with B-esterase type was reported from diamondback moth (Mohan and Gujar, 2003). In addition to qualitative changes in esterase banding pattern, there is a significant increase in the production of the CEs in at least three field populations. The over-production and appearance of new esterase isozymes might be the consequence of increased transcription or gene duplication events (Malathi et al., 2015;Panini et al., 2016).
Most of the GST gene expansions are from epsilon and delta subclasses together accounting for 18 genes (54.5%) out of 33 genes identified from the draft genome of L. orbonalis. These two GST subclasses are insect-specific, widely expanded in many other insect species whereas other classes have wider taxonomic distribution (Friedman, 2011). There were at least nine GSTs over-expressed >50-fold in one or more of the field-collected resistant populations of L. orbonalis. In Anopheles gambiae, Apis mellifera, B. mori, D. melanogaster, Nasonia vitripennis, and T. castaneum respectively 28,8,23,37,19,35 GSTs and 51,24,76,35,41,49 CEs were identified (Yu et al., 2008(Yu et al., , 2009Oakeshott et al., 2010). Similarly, the gene expansion of cytochromoe p450 monooxygenases (CYP) and their role in insecticide resistance in L. orbonalis were recently documented (Kariyanna, 2019;Kariyanna et al., 2020a). Together with the results of present and earlier studies it is very much evident that multiple metabolic enzymes contribute toward overall insecticide resistance in L. orbonalis. However, in many insect species including L. orbonalis, it is difficult to distinguish whether the expansion of these metabolic genes and their involvement in insecticide detoxification are driven from ancestry or adaptation to stress or both.

CONCLUSION
In conclusion, the present study provides the first genomic level information on GST and CE gene families of L. orbonalis and their role in insecticide resistance. The information further supports insect resistance management through use of reduced risk pesticides including biologicals.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
BKr: conducted the experiment, analysis, manuscript writing, materials and methods, and hypothesis. AP: analysis, material methods, and corrections. RA: analysis and technical suggestion, manuscript correction, and resource supply. AA and PJ: bioinformatics and analysis. RG and BKl: resources supply and, technical guidance and suggestions. TV: resources supply, technical guidance and suggestions, and reviews and literatures. MB and JD: analysis, technical guidance and suggestions, and reviews and literatures. MM: designing experiment, resources supply, manuscript corrections, analysis, material and methods, technical guiding, and suggestions. YP: technical guidance and suggestions, and reviews and literatures. All authors contributed to the article and approved the submitted version. necessary support. Authors are thankful to Director, ICAR-NBFGR for timely release of fund. The facilities extended by Director, ICAR-NBAIR, and Director, ICAR-IIHR, Bengaluru are gratefully acknowledged.