Functional analysis of CYP6ER1, a P450 gene associated with imidacloprid resistance in Nilaparvata lugens

The cytochrome P450 CYP6ER1 has been reported to play an important role in imidacloprid resistance of the brown planthopper (BPH), Nilaparvata lugens, and is overexpressed in most resistant populations. In the present study, we confirmed that CYP6ER1 expression can be induced by certain levels of imidacloprid. Developmental expression analysis revealed that CYP6ER1 was expressed highly in the adult stage, and tissue distribution analysis showed that CYP6ER1 was expressed mainly in the fat body and midgut. RNA interference (RNAi) of CYP6ER1 and transgenic expression of CYP6ER1 in Drosophila melanogaster both suggested that the expression of CYP6ER1 is sufficient to confer imidacloprid resistance. Furthermore, we analyzed the interaction of imidacloprid and CYP6ER1 monooxygenase by using dynamic simulations and molecular docking. We found that Nitrogen atoms in the heterocycle of the imidacloprid molecule may bind to iron atoms in the center of the homology model of CYP6ER1 via 4,5-dihedro-1H-imidazole. This finding contributes to a better understanding of how CYP6ER1 takes part in the insecticide metabolism.

Scientific RepoRts | 6:34992 | DOI: 10.1038/srep34992 structural data when experimental techniques fail. Many proteins are too large for NMR analysis and cannot be crystallized for X-ray diffraction 22 , resulting in a lack of elucidate crystal structures for P450s and their substrate complexes. In these cases, homology modeling gives an insight into understanding the complexity of insecticide resistance mediated by P450 genes 23,24 . By comparing homologous candidates, a number of candidates searched in the database can, to a certain extent, allow researchers to predict and visualize specific molecular structures and interactions, which are difficult to prove using traditional biological experiments. Jones utilized homology modeling based on a template of the X-ray structure of the phylogenetically related human CYP3A4, to simulate the interaction between D. melanogaster CYP6G1 and DDT and neonicotinoids, and indicated roles for CYP6G1, CYP12D1 and CYP6A2 in flies overexpressing P450 enzymes with broad substrate specificities in chemical metabolism 25 . Karunker found interactions of imidacloprid with the biotype Q variant of the CYP6CM1 enzyme might involve key amino acids such as Phe-130 and Phe-226 in the enzyme active site with the lowest energy, and they predicted a putative hydroxylation site 26 . Unfortunately, while the application of homology modeling is now widely spread in the medical sciences, it has been used in only a limited capacity for the determination of insecticide metabolism mechanisms in pest control.
ln this study, we determined, using RNAi and transgenic expression in Drosophila, that the expression of P450 CYP6ER1 is sufficient to confer imidacloprid resistance. Additionally, we utilized molecular modeling in order to predict the catalytic site geometrics and identify the key residues responsible for imidacloprid binding with CYP6ER1. These experiments contribute to clarifying the enzyme responsible for metabolizing imidacloprid and causing the resistance phenotype.

Results
CYP6ER1 expression patterns in response to insecticide exposure. In order to confirm the expression pattern of CYP6ER1 in insecticide resistant BPH strains, one-day-old female adults of a field-caught strain from Guangxi and a susceptible contrasted strain from Zhejiang were selected and the field-caught strain was further treated with 4 different doses of imidachloprid. We found that the mRNA level of CYP6ER1 in the fieldcaught strain was 2.1-fold higher than that observed for the susceptible strain, and the mRNA expression level in the samples treated with 30 ppm and 60 ppm imidachloprid was significantly higher than that observed for the untreated field-caught strain. However, there was no significant difference between the untreated fieldcaught strain and samples treated with 90 ppm or 120 ppm imidachloprid (Fig. 1). These results suggested that the up-regulation of CYP6ER1 can be triggered by imidachloprid, and it can reach at maximum within certain dosages of imidachloprid.
Developmental expression and tissue distribution of CYP6ER1. Quantitative Real-time PCR (qRT-PCR) was used to detect the expression patterns of CYP6ER1 during N.lugens from the first instar nymphs to the 13th day adulthood. The mRNA levels developed a stable trend of low expression before 5th instar, but increased significantly from adult day 1 to a relatively high level between adult day 9 and day 13 ( Fig. 2A). We also sampled different tissues from 1-day-old female adults to investigate their mRNA levels of CYP6ER1, and it was expressed highly in the fat body and midgut, expressed at a relatively lower level in the ovaries, Malpighian tubule Figure 1. CYP6ER1 expression patterns in response to insecticide exposure. Quantitative real-time PCR measurement of the fold change in the expression of CYP6ER1 in the field-caught strain compared with the susceptible strain, while the field-caught strain treated with 4 different doses of imidacloprid (30ppm,60ppm,90ppm and 120ppm, respectively) compared with the field-caught strain without imidacloprid. β-actin was used as an internal reference gene and the mRNA expression level of CYP6ER1 in the susceptible strains was designated as one. The data represent the mean + SEM (n = 30). *P < 0.05 level (t-test) and **P < 0.01 level (t-test).
Scientific RepoRts | 6:34992 | DOI: 10.1038/srep34992 and cuticle, and not expressed in the muscle (Fig. 2B). The high levels of expression in the essential metabolic organs correspond to the physiological role of CYP6ER1 in imidacloprid resistance.

RNA interference of CYP6ER1.
To evaluate the contributions of CYP6ER1 in imidacloprid resistance in vivo, we designed dsRNA to silence CYP6ER1 in female adult from the Guangxi strain. The mRNA expression of CYP6ER1 decreased within 48 hours after injecting three doses (50 ng, 100 ng and 200 ng) of dsCYP6ER1, but there was no significant difference between the three doses of dsRNA (Fig. 3A). We chose 50 ng imidacloprid treatment for 48 hours to test BPH mortality after RNAi (Fig. 3B). The mortality of the dsCYP6ER1 sample was 40.00% and 34.42% higher than those of the ddH 2 O and dsGFP samples, respectively. This suggests that the expression of CYP6ER1 should be of great importance to imidacloprid resistance in BPH.
Transgenic expression of the potential resistance genes CYP6ER1 in D. melanogaster. To identify whether the expression of CYP6ER1 is sufficient to confer imidacloprid resistance, we used a transgenic approach utilizing the GAL4/UAS system of D. melanogaster. CYP6ER1 was inserted into the vector pUAST and embryonic microinjections were performed on D. melanogaster w 1118 to generate the transgenic germline carrying a UAS-CYP6ER1 transgene. The flies subsequently were crossed to the Act5C-GAL4 line. The mRNA levels of CYP6ER1 were determined in the progeny (1-3 days-old adults) by RT-PCR, confirming the expression of the transgene in the Act5C > CYP6ER1 line (progeny expressing the CYP6ER1) (Fig. 4A). In comparison, the CYP6ER1 mRNA was not detected in another two control D. melanogaster lines (UAS-CYP6ER1 and Act5C-GAL4, progeny not expressing the CYP6ER1).
Bioassays showed that Act5C > CYP6ER1 line was resistant to imidacloprid at a treated dose of 30 μ g/ml imidacloprid, with significant higher survival rate (86.42%, p < 1e-10 Chi-squared Test) than the control lines (4.23% and 11.15% respectively) (Fig. 4B). These results suggest that the expression of CYP6ER1confers imidacloprid resistance.
Homology modelling of CYP6ER1. Structure based sequence alignment was used to compare the key residues previously predicted to be responsible for imidacloprid binding and metabolism, such as CYP6CMA1,  Table 1). CYP6ER1 had a high similarity and shared several functionally identical amino acids, suggesting that CYP6ER1 may have the same functional domain as other imidacloprid metabolism proteins. The prediction of the substrate recognition site (SRS) secondary structure elements and key residues were cited in the research of CYP6CM1 26 . A BLAST search with the default parameters for CYP6ER1 in the PDB_ nr95 database identified the homologous protein, CYP3A4 (PDB code: 3NXU), with an intermediate sequence identity of 31%, which was then used to construct the CYP6ER1 structure.
A total of 20 models were generated. Among these 20 models, 6ER1.M0011 had the smallest value of PDF (probability density function) Total Energy and PDF Physical Energy, and relatively low DOPE Score (Table S3), by using Cut overhangs, refine loops and use of DOPE method as parameters. This indicated that 6ER1.M0011 The mRNA expression level of CYP6ER1 in dsGFP ground and dsCYP6ER1 group were relative to β -actin mRNA level and the mRNA expression level of CYP6ER1 in dsGFP group was designated as one. (B) Effects on the mortality rate after injection of water or dsGFP or dsCYP6ER1 and followed by application of imidacloprid (50 ng/ ml). Double asterisks indicate significant differences in the mRNA relative expression level between the treatment and control groups (P < 0.01, t-test). The data represent the mean + SEM (n = 3), **P < 0.01 level (t-test). The comparison between survival rates of two control lines and transgenic line exposed to 30 μ g/ml imidacloprid. The data shown are the mean ± SEM (n = 3). **P < 0.01 (Chi-squared Test).
was the most suitable model among the 20 counterparts, and could be utilized to simulate the molecular docking. 6ER1.M0011 was then analyzed by Ramachandran and most of these residues were located in the most favored region, with 11 out of 470 residues in the disallowed region, suggesting it was a relatively good quality model but needed further regulation and optimization.
The binding mode of imidacloprid to the CYP6ER1 active site. During our analysis, we found ritonavir as specific inhibitor of CYP3A4 in human, which suggested that its docking method may be similar to   Table 1. Alignment of the analogous imidacloprid binding residues in four imidacloprid metabolizing enzymes: CYP6ER1, CYP6CM1, CYP3A4, CYP6G1. Analyses were based on structure sequence alignment between the CYP6CM1 model and the CYP3A4 crystallographic structure together with sequence alignment to CYP6G1. Some of the residues that contribute to the stabilization of imidacloprid in the active site are also important in the stabilization of the heme and oxygen molecules.
CYP6ER1. Imidachloprid and ritonavir both have a similar heterocycle (Fig. S1A,B) and ESP charge distribution (Fig. S1C,D). Rritonavir bonds covalently to CYP3A4 via nitrogen atoms, and this mode may be similar to the way imidachloprid bonds to CYP6ER1.
We used the homology model based on CYP6ER1 as the receptor, and the imidachloprid molecule as the ligand. Nitrogen atoms in the heterocycle of the imidachloprid molecule bonded to iron atoms in the center of the homology model via 4,5-dihedro-1H-imidazole. The temperature of CYP6ER1-imidacloprid compound hovered around 300 K, with its deviation below 10 K (Fig. S2), and the total energy, kinetic energy and potential energy started from high energy to a stable relatively low level after about 5 picoseconds (Fig. S3). RMSD of atomic coordination of the protein backbone, HEME and ligand fluctuated up and down in the range between 1.20 Å to 1.50 Å (Fig. S4). The results confirmed the hypothesis of the formation of a CYP6ER1-imidacloprid complex based on its relatively credible stability and balance of its energy and the structure. Finally, the interaction energy of CYP6ER1-imidacloprid complex was calculated at − 114.66 kcal/mol, indicating that this may be the ideal model to simulate the P450-HEME-IM4 complex and its binding pocket, and suggested that some amino acids in the center of CYP6ER1 have a role in the molecular docking the between CYP6ER1 and imidachloprid (Figs 5 and 6 and Table 1). Most of these residues were predicted to generate the hydrophobic interface of the binding cavity, in which Phe-124 was predicted to anchor imidacloprid by aromatic interaction generating the largest hydrophobic patch in the interface. We predicted that imidacloprid would be further stabilized by hydrogen donors and acceptors that include Pro-375 (SRS-5) and Pro-289 (SRS-2), which can interact with imidacloprid's imine moiety and Leu-312(SRS-4).

Discussion
It is generally believed that detoxifying enzymes, referred to as esterase, glutathione S-transferase and P450s, contribute to defending against and catabolizing toxins and insecticides in insects, and that many involve point mutations of target genes for insecticide which increase the copy number of the gene, mRNA level and the diversification of the coding sequence 27 . In a study of neonicotinoids, total enzyme activities of the candidate metabolic enzymes were measured in N. lugens, indicating that the increase of neonicotinoids resistance was due mainly to P450 monooxygenase, instead of mutations of the nicotinic acetylcholine receptors that were considered as possible targets 28 . Puinean used piperonyl butoxide, an inhibitor specific to P450s, to decrease the level of imidacloprid resistance significantly, which confirmed the important roles of P450s in insecticide resistance 29 . Similarly, the phenomenon was also reported in aphids 30 . Additionally, specific P450 genes have been related to neonicotinoids resistance. For instance, Drosophila CYP6G1 has been proven to confer imidacloprid resistance 14,31 , and Puinean suggested that CYP6Y3 can confer neonicotinoids resistance to Myzus persicae 32 .
The P450 CYP6ER1 has been associated with imidacloprid resistance in N. lugens 19 . Here, we found that the mRNA level of CYP6ER1 was significantly higher in the field-caught resistant strains than in the susceptible strains, which generally coincides with previous references. Susceptibility to imidacloprid increased by 75.61% after RNAi treatment in the experimental population (Fig. 3). Similarly, it has been reported that susceptibility to Deltamethrin in Locusta migratoria increased by 92.14% and 72.93% after RNAi knockdown of CYP409A1 and CYP409B1, respectively. Additionally, transgenic expression of CYP6ER1 in D. melanogaster lead to a significant increase of imidacloprid resistance. We therefore conclude that the expression of CYP6ER1 is sufficient to confer imidacloprid resistance. Of particular note is that fact that, though the expression of CYP6ER1 in the resistant strain can be induced by imidacloprid treatment, the over-expression of this gene was limited when BPH were treated by a relatively higher dose of imidacloprid. This indicates that the imidacloprid resistance phenotype may Figure 6. Overview of P450-HEME-IM4 complex and its binding pocket. The one conformation of MDbased optimized structure. Imidacloprid is indicated with light blue carbon atoms and the HEME is showed with grass green atoms. Predicted binding residues within 4 Å of imidacloprid are showed in pink. Phe-124 and Asp-290 anchor imidacloprid by aromatic interactions, and generate significant part of the contact surface area. Imidacloprid may be further stabilized by hydrogen bonded with Pro-375 and Pro-289.
Scientific RepoRts | 6:34992 | DOI: 10.1038/srep34992 be more complex than just the over-expression of a single gene. Thus, it is necessary to review and optimize the traditional strategy of pest control.
P450 monooxygenase is well known for metabolizing xenobiotics, such as insecticides, to less toxic metabolites 33,34 . For example, CYP6A1, 6A2, 6A5, and 12A1 all have some manner of cyclodiene epoxidase biochemical activity. Graham and Peterson observed the general topology and tridimensional fold of all cytochromes P450 35 . Only a few strictly conserved residues, as components of a four-helix (D, E, I and L) bundle, two sets of β -sheets, and a coil, participate to form the center of P450s together with a central heme-group bound to the thiolate of a relatively highly conserved cysteine residue, which give access to the biotransfer of electrons and protons to activate oxygen. On the other hand, the individual substrate specificity of the enzymes lies on the individual residues despite their homologous and similar regions, conferring their unique function in reactions with different insecticides.
Over-expression of CYP6ER1 was able to confer imidacloprid resistance, and this could be induced by certain doses of imidacloprid itself. Therefore, it is important to explore further the interaction between CYP6ER1 and imidacloprid. Utilizing homology modeling technology, and based on the high degree of homology with CYP3A4 in human, we predicted and visualized the molecular interaction wherein nitrogen atoms in the heterocycle of the imidachloprid molecule bond mainly to iron atoms in the center of the homology model via 4, 5-dihedro-1H-imidazole. CYP6CMAvQ and CYP6G1, two identified P450s which detoxify imidachloprid, were referred to our alignment with CYP6ER1 and CYP3A4 26,36 , confirming that target sequences and the substrate recognition site (SRS) regions had a high level of similarity 26 . Alignment of analogous imidacloprid binding residues in these four imidacloprid metabolizing enzymes indicated a number of consistent binding residues among them. To assess the stability and the balance of the CYP6ER1-imidacloprid complex, several parameters in the field of energetic and geometric criteria were calculated and analyzed until we obtained an optimized model to simulate the interaction. After optimization, our homology model proved to be of high quality. Additionally, there were some putative key residues responsible for imidacloprid binding with similarity among the four enzymes, most of which create the hydrophobic interface of the binding cavity, as suggested by the docking analysis. Phe-124 and Asp-290 anchor imidacloprid by aromatic interactions, and create a significant part of the contact surface area between imidacloprid and CYP6ER1. Pro-375 (SRS-5) and Pro-289 (SRS-2) likely generate hydrogen bonds to further stabilize imidacloprid. In these residues interacting with imidacloprid, Ala-313 and Gly-314 are conserved in three other imidacloprid metabolizing enzymes: CYP6CM1, CYP3A4 and CYP6G1. The conservation of these residues and the relative large hydrophobic contact surface area, show that residues play an important role in contributing to these enzymes function of metabolizing imidacloprid. However, we cannot ignore the possibility of other residues constructing the model, despite the lack of apparent homologues among these enzymes, which may develop its own specific function and structure for the enzyme.
Our findings have demonstrated the importance of homology modelling and docking to search the potential catalytic center of an enzyme with its substrate. Our research can also lead a strong foundation for future insecticide design. However, a deeper understanding of the interaction mode of imidacloprid related heterocyclic compounds with P450s is needed in order to guide the rational development of insecticide-related compounds and increasing the accuracy and efficiency of computational based drug construction. In addition, the predicted residues and sites involved in detoxifying imidacloprid require further biochemical validation.

Material and Methods
Insect strains. The field-caught strains were collected from Nanning, Guangxi, China in 2012, and provided by Institute of Plant Protection, Guangxi Academy of Agricultural Sciences. The susceptible strains were obtained from Zhejiang Academy of Agricultural Sciences in 2012 and reared on rice plants without insecticide treatment. The resistance ratios for imidacloprid of the field-caught strain (LC 50 249.13 mg/L) to the susceptible strain (LC 50 0.22 mg/L) was 1132-fold.
Insecticide exposure. Imidacloprid (96.04%) was purchased from Tianyuan Biochemical Co. Ltd (Guangxi, China). We use the rice stem dipping method for insecticide exposure to BPH as described 37 . Different doses of imidacloprid were diluted with acetone to a series of concentrations, and then 1% Tween 80 was added. Rice stems rooted in culture cups were dipped in insecticide solutions for 30 s. Once the stems were dry, one-day-old long-winged female adults were put into each culture cup and fed on these treated rice. Each treatment was performed in triplicate.
Total RNA extraction. Total RNA was extracted from BPH by using a Total RNA extract kit (Omega).
Genomic DNA was digested with DNase I (Takara Bio Inc., Kyoto, Japan). The quality and concentration of RNA samples were examined by agarose gel electrophoresis and spectrophotometric analysis. cDNA was synthesized by reverse transcription in 20 μ l reaction volumes containing 1 μ g of total RNA by using PrimeScript reverse transcriptase (Takara). The reaction mixture was stored at − 20 °C for future use 37 . qRT-PCR. Primers were designed for qRT-PCR based on the complete coding sequences of CYP6ER1 genes as archived in the NCBI database (Table S1). The primers had annealing temperatures of approximately 58 °C. Primers corresponding to the housekeeping gene β-Actin were used as endogenous controls 38 Table S1. The dsRNA synthesis step was performed as described previously 37 . For the RNAi efficiency evaluation, 0.05, 0.1 and 0.2 micrograms of dsCYP6ER1 were injected into the side of the abdomen of synchronous one-day-old female N.lugens macropterous using a 10 μ l micro-syringe (Hamilton). The control injection was performed with an equivalent volume of dsGFP. Each experimental group had 30 individuals. Ten adults were selected randomly at 24 h, 36 h and 48 h after injection for the independent detection of CYP6ER1.
Total RNA extraction and qRT-PCR detection were performed as described above. For the susceptibility detection after RNAi, thirty adults from each group were chosen 24 h after injection, with 0.05 micrograms of ddH 2 O, dsGFP and dsCYP6ER1 for each one day-old macropterous female adult. Thirty adults were placed in rice soaked in a solution of 240 ng/μ l imidacloprid for 30 seconds. The N.lugens mortality was calculated after 96 hours. Each experiment was repeated in triplicate.
Developmental expression and tissue distribution analysis. The expression of CYP6ER1 were estimated by qRT-PCR as described 38 from the first instar nymphs to 13-day-old female adults from the field-caught strain without imidacloprid. To investigate the expression of CYP6ER1 in different tissues, we extracted the total RNA from the fat body, ovary, midgut, muscle, Malpighian tubule and cuticle. The following steps were performed as described 38 .
Construction of Transgenic Drosophila and Bioassays. The CYP6ER1 CDS cloned was inserted into a pValium20 vector to prepare UAS-CYP6ER1 constructs. Subsequently, the constructs were microinjected into the embryos of the y sc v nanos-integrase; attP2 D. melanogaster germline using standard techniques in the Center of Biomedical Analysis, Tsinghua University. The transformed line was crossed with anAct5C-GAL4 line for expression of the CYP6ER1 gene. The genotype of the cross was Act5C > CYP6ER1. For use as a control, the transformed line and the Act5C-GAL4 line were crossed with the W 1118 line, the progenies of which did not express CYP6ER1 gene with the genotype of UAS-CYP6ER1 and Act5C-GAL4. RT-PCR was used to confirm the expression of the CYP6ER1 gene in transgenic Drosophila using primers specific for the CYP6ER1 gene and the reference housekeeping gene DmActin 40 . For each cross, total RNA was extracted from 10 adults using the methods described above and three biological replicates were performed for each experiment. The reaction conditions were 94 °C for 5 min, 30 cycles of 94 °C for 30 s, 58 °C for 30 s, and 72 °C for 30 s, using GoTaq G2 Green Master Mix (Promega). Twenty 1-3 day old adults (10 females and 10 males) were used in the contact bioassays, as described previously 41 . In brief, flies were placed in vials with 10 ml corn meal medium containing 300 μ g imidacloprid. Three replicates were performed for each assay. After 72 h, the survival rates were calculated and analyzed using the Chi-squared test.
Homology modelling. The structure of CYP6ER1 was predicted using the Homology Modeling protocol of Discovery Studio 2.5. A BLAST search against the PDB_nr95 database was used to search the homology sequence by using the E-value 10-10 as a parameter. A multi sequence alignment between CYP6ER1 and homology sequences was performed to obtain the best reference for the homology model building. An ideal homologous model should cover as much of the Alignment Length of our target sequence as possible, while keeping relatively high Positive and Sequence identity, which was also analyzed via Align Multiple Sequence in Discovery Studio with default parameters. The crystal structure of human CYP3A4 with the highest resolution and good stereochemistry was used to build the homology model. The best model for subsequent optimization was chosen according to the total PDF energy, PDF physical energy and DOPE score, using Cut overhangs, refine loops and the use DOPE method as parameters.

Molecular dynamics simulations.
According to the covalent bonds between ritonavir (RIT) and CYP3A4, we speculated that bonding also occurred between imidacloprid (IM4) and CYP6ER1, though conventional docking programs cannot be utilized to generate the ligand-bound conformations. Observing charge distribution in local structures of imidacloprid and CYP6ER1 calculated by the Gaussian-09 program at the HF/6-31G* level (Gaussian 09, revision a.02. Gaussian Inc., Wallingford) 42 , we found 4, 5-dihedro-1H-imidazole in structural imidacloprid was partly similar to the thiazole structures of RIT. Therefore, it is possible that nitrogen atoms in the heterocycle of the imidachloprid molecule bonded to iron atoms via 4, 5-dihedro-1H-imidazole. A comparison of the electrostatic distributions of RIT and IM4 suggests that IM4 likely has the same binding mode as RIT, and IM4 was found to be manually docked into CYPYER1. A molecular dynamics (MD)-based optimization assessment of the complex structure of CYPYER1 protein was performed with the Discovery Studio 2.5 Dynamics (Heating or Cooling, Equilibration and Production) Module under the CHARMm force field 43 . The partial charge and implicit solvent model were set to Momany-Rone and generalized born as the molecular dynamics (MD)-based optimization parameters, Electrostatics as Spherical cutoff, Apply SHAKE constraint as True, Minimization as 400, Heating as 2,000, Equilibration as 2,000, Production as 100,000 and other parameters as the defaults.