Genomic and transcriptomic insights into the cytochrome P450 monooxygenase gene repertoire in the rice pest brown planthopper, Nilaparvata lugens

YGENO-08762; No. of pages: 9; 4C: Genomics xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Genomics journal homepage: www.elsevier.com/locate/ygeno

Genomic and transcriptomic insights into the cytochrome P450 monooxygenase gene repertoire in the rice pest brown planthopper, Nilaparvata lugens Shu-Hua Lao 1, Xiao-Hui Huang 1, Hai-Jian Huang, Cheng-Wen Liu, Chuan-Xi Zhang, Yan-Yuan Bao ⁎ State Key Laboratory of Rice Biology and Ministry of Agriculture Key Laboratory of Agricultural Entomology, Institute of Insect Sciences, Zhejiang University, Hangzhou 310058, China

a r t i c l e

i n f o

Article history: Received 10 April 2015 Received in revised form 26 June 2015 Accepted 9 July 2015 Available online xxxx Keywords: Planthopper Cytochrome P450 monooxygenase Insecticides

a b s t r a c t The cytochrome P450 monooxygenase (P450) gene family is one of the most abundant eukaryotic gene families that encode detoxification enzymes. In this study, we identified an abundance of P450 gene repertoire through genome- and transcriptome-wide analysis in the brown planthopper (Nilaparvata lugens), the most destructive rice pest in Asia. Detailed gene information including the exon–intron organization, size, transcription orientation and distribution in the genome revealed that many P450 loci were closely situated on the same scaffold, indicating frequent occurrence of gene duplications. Insecticide-response expression profiling revealed that imidacloprid significantly increased NlCYP6CS1v2, NLCYP4CE1v2, NlCYP4DE1, NlCYP417A1v2 and NlCYP439A1 expression; while triazophos and deltamethrin notably enhanced NlCYP303A1 expression. Expression analysis at the developmental stage showed the egg-, nymph-, male- and female-specific expression patterns of N. lugens P450 genes. These novel findings will be helpful for clarifying the P450 functions in physiological processes including development, reproduction and insecticide resistance in this insect species. © 2015 Elsevier Inc. All rights reserved.

1. Introduction The brown planthopper, Nilaparvata lugens Stål (Hemiptera: Delphacidae) is the most devastating rice plant pest throughout Asia. N. lugens causes extensive damage to rice by sucking rice phloem sap and blocking phloem vessels, which stunts plant growth [1]. As with many other sucking pests, N. lugens is an efficient vector of plant viruses, including rice ragged stunt virus and rice grassy stunt virus, which can cause rice ‘grassy stunt’ and ‘ragged stunt’ diseases, respectively [2]. Together, these factors can cause severe yield reductions and significant economic losses [3]. Currently, chemical control using a variety of insecticides, including organophosphates, carbamates, pyrethroids and neonicotinoids, remains the first choice for managing N. lugens populations in Asia; however, the widespread use of insecticides has given rise to planthopper resurgence and additional environmental risks [4]. The neonicotinoids, which are highly efficient and exert long-lasting effects, have become the most widely used insecticides for N. lugens control since the introduction of the first commercialized imidacloprid by Bayer in the early 1990s [3]. However, N. lugens field populations have developed a high level of resistance to imidacloprid over large areas in China since 2005, which has led to control failure [5,6]. In N. lugens, imidacloprid ⁎ Corresponding author. E-mail address: [email protected] (Y.-Y. Bao). 1 These authors contributed equally.

resistance was attributed to enhanced oxidative detoxification by over-expressed P450 monooxygenases (P450s) [7–9]. Two N. lugens P450s, CYP6ER1 and CYP6AY1, have been elucidated to be involved in imidacloprid resistance [3,10]. The discovery of functional P450s is now an area of intense interest because of potential for these candidate molecules to mediate insecticide resistance. P450s form the largest and most diverse superfamily of multifunctional enzymes that play important roles in the biosynthesis of endogenous compounds and detoxification of xenobiotic compounds [11]. Insect P450s can be classified into four major clades: three microsomal cytochrome P450 clades (CYP2–CYP4) and the mitochondrial cytochrome P450 clade. In insects, P450s contribute to numerous functions. They are implicated in biosynthetic pathways of juvenile hormones and ecdysteroids, which are at the center stage of insect growth, development and reproduction [12–16]. Some P450s are involved in metabolism of plant secondary metabolites, which is important for the adaptation of insect herbivores to their plant hosts. Some P450s are responsible for the oxidative detoxification of synthetic insecticides, resulting in metabolic resistance to insecticides [17,18]. Recently, we completed the whole-genome sequencing of N. lugens and obtained the gene annotation (Zhejiang University N. lugens genome project team) [19]. The N. lugens genome is the first characterized genome of a monophagous sap-sucking arthropod herbivore, and the availability of its genome sequence allows the global analysis of P450 genes, which will yield insights into the resistance mechanisms against a variety of insecticides in this hemimetabolous insect species. Previously, we

http://dx.doi.org/10.1016/j.ygeno.2015.07.010 0888-7543/© 2015 Elsevier Inc. All rights reserved.

Please cite this article as: S.-H. Lao, et al., Genomics (2015), http://dx.doi.org/10.1016/j.ygeno.2015.07.010

2

S.-H. Lao et al. / Genomics xxx (2015) xxx–xxx

performed transcriptome sequencing using next-generation highthroughput Illumina technology, which generated detailed transcript information about the developmental stages, wing dimorphism, sex differences, immune responses and tissue specificity in N. lugens [20–22]. Here, a thorough search of the N. lugens genome sequence coupled with analyses of the available transcriptome datasets provided comprehensive information about the P450 genes, which can be classified into at least four P450 subgroups based on sequence homology and subcellular location. A high-throughput RNA sequencing (RNA-Seq) analysis revealed the expression specificity of P450 genes in response to different insecticide treatments in N. lugens. These data might be helpful for clarifying the functional roles of N. lugens P450s and in establishing their association with insect resistance to chemical insecticides. 2. Results and discussion 2.1. Identification of the N. lugens P450 genes in the genomic and transcriptomic databases We identified 68 predicted P450 gene sequences by searching the N. lugens genomic and transcriptomic databases. These sequences were confirmed using the tBLASTX algorithm with a cut-off E-value of 10−10 (Table 1). We retrieved 59 sequences that best matched with sequences of P450 homologs of the small planthopper Laodelphax striatella. We identified nine sequences with high similarity to the P450s of several insect species, including Sogatella furcifera, Acyrthosiphon pisum and Aedes aegypti. The N. lugens P450 genes contain complete open reading frames that each encodes ~ 500 amino acids (Table 1). We submitted all N. lugens P450 gene sequences to David Nelson (http://drnelson.uthsc.edu/CytochromeP450.html) in the P450 nomenclature committee and obtained official nomenclature of these P450s. The N. lugens P450 genes can be classified into four major clades. The CYP4 clade is the largest, and is followed by the CYP3 clade. The N. lugens CYP4 clade includes 14 NlCYP4, one NlCYP380, five NlCYP417, three NlCYP425, one NlCYP426 and three NlCYP439 genes. The members of this clade in some insect species have been linked to xenobiotic metabolism and insecticide resistance. The high diversity of N. lugens P450 genes in the CYP4 clade is reflected in a great diversity of gene functions. The N. lugens CYP3 clade is also a large subfamily that includes 13 NlCYP6, one NlCYP408, one NlCYP418, three NlCYP427 and one NlCYP3115 gene. However, the genes for CYP6a2, CYP6a13, CYP6a14, CYP6k1 and CYP3A18 in the CYP3 clade were all absent in N. lugens. Some insect P450s in this clade appear to be inducible metabolizers of xenobiotics, while others have been linked to odorant or pheromone metabolism. The N. lugens gene losses in this clade might be associated with its reduced need for detoxification or recognition needs for combatting its sole plant host, rice. In the CYP2 clade, 10 genes, including NlCYP15, NlCYP18 and NlCYP303–307, were identified in the N. lugens genome. CYP306 and CYP307 genes have been found to be involved in insect ecdysone metabolism during developmental stages [16,23]. In some insect species, CYP15 is the stereospecific epoxidase involved in juvenile hormone biosynthesis [24], while CYP18 is known to be an ecdysone-responsive gene, whose expression was increased by the molting hormone [12]. We compared the currently available insect species genome sequences and found that CYP18 is absent in the mosquito Anopheles gambiae and CYP15 is absent in the fly Drosophila melanogaster, but these two genes were present in N. lugens. The mitochondrial P450s are not found in plants or fungi, and appear to be restricted to animals. The N. lugens mitochondrial P450 family includes NlCYP301, NlCYP302, NlCYP314, NlCYP315, NlCYP353, NlCYP404 and NlCYP419 genes. Among the mitochondrial P450 genes, CYP404 homologs in insect species have only been identified in the small planthopper L. striatella and a termite Zootermopsis nevadensis. In vertebrates, mitochondrial P450s are all involved in steroid or vitamin D metabolism. In insects, only a few mitochondrial P450s have annotated physiological functions. For example, a CYP12

gene was first identified in the housefly because of its importance in the metabolism of a variety of xenobiotics. The Drosophila and Anopheles CYP12 genes have been associated with insecticide resistance; however, N. lugens lacks a CYP12 gene. The reduction in P450 genes in N. lugens seems to be correlated with the levels of exposure to external stressors, such as host defense systems, environmental toxins and a single plant host. Analysis of exon–intron organization revealed that most N. lugens P450 genes consist of multiple exons (Table 1). Some genes were located on the same scaffolds and formed gene clusters. For example, the following genes were on the indicated scaffolds in the same transcriptional orientation: NlCYP6AX1v2, scaffold297; NlCYP6CS1v2 and NlCYP6FK2, scaffold53; NlCYP418A1v2 and NlCYP427A1, scaffold3261; NlCYP427A1 genes, scaffold59; NlCYP417B1 genes, scaffold926; and NlCYP4C genes (including NlCYP4C61v2, NlCYP4C76 and NlCYP4C77), scaffold 745 (Fig. 1 and Table 1). Additionally, the following genes were found on the same scaffold but in opposite orientations: NlCYP18A1 and NlCYP306A2 genes, scaffold953; and NlCYP417 genes (including NlCYP417A1v2, NlCYP417A2 and NlCYP417A3), scaffold371 (Fig. 1 and Table 1). The fact that two or more P450 loci were located on the same scaffold implies that N. lugens has undergone many gene duplication events in its genome. Among the hemimetabolous sapsucking insect species, only the A. pisum and N. lugens genomes have been sequenced to date. In the A. pisum genome, several assembled scaffolds contain two or more P450 loci that are closely related to each other [25] — this observation was thought to be specific to the A. pisum genome [18]. Here, our findings provide new evidence for P450 gene evolution through tandem duplication events in N. lugens. Phylogenetic analysis showed that N. lugens P450s evolved into four major clades: CYP2–CYP4 and mitochondrial P450 (Fig. 2). The genes of the CYP3 clade were more closely related to those in the CYP2 and CYP4 clades than to those in the mitochondrial P450. The N. lugens P450s showed the closest phylogenetic relationship with homologs of two hemimetabolous species, Laodelphax striatellus and S. furcifera. Among the N. lugens P450 genes, NlCYP4CE1v2 showed significant sequence similarities to several hydrobios, including Daphnia pulex, Daphnia magna, Carcinus maenas and Portunus trituberculatus, but low similarities to other insect species. This P450 gene may have unique functions in N. lugens. 2.2. Construction of RNA-Seq cDNA libraries and confirmation of the insecticide-responsive P450 genes To better characterize gene expression changes in response to different insecticide treatments, we constructed RNA-Seq cDNA libraries generated from imidacloprid-, triazophos- and deltamethrin-treated 5th instar nymphs using the Illumina HiSeq™ 2000 sequencing platform. Acetone-treated specimens were used as a control. We performed a classification analysis aimed at identifying differentially expressed genes with more than twofold expression changes between insecticide- and acetone-treated samples by comparison with the Gene Ontology (GO) annotation database (http://wego.genomics.org.cn/cgi-bin/wego/index. pl). The N. lugens genes, which showed significant changes in expression levels in response to imidacloprid treatment, were categorized into 43 functional groups (Fig. 3A), while genes with the most significant expression changes in response to triazophos and deltamethrin treatments were categorized into 24 and 18 functional groups, respectively (Fig. 3B). In each of the three main categories (biological process, cellular component and molecular function) of the GO classification scheme, the ‘Metabolic process,’ ‘Cell or Cell part’ and ‘Catalytic activity’ terms were the most commonly associated with the three comparisons, implying that the expression levels of these genes were most affected by the three insecticides. We compared the differentially expressed genes involved in metabolic processes. The criteria of false discovery rate (FDR) ≤ 0.001 and log2 ratio ≥ 1 were used to identify the most significantly differentially expressed genes between the insecticide- and acetone-treated

Please cite this article as: S.-H. Lao, et al., Genomics (2015), http://dx.doi.org/10.1016/j.ygeno.2015.07.010

S.-H. Lao et al. / Genomics xxx (2015) xxx–xxx

3

Table 1 Identification of P450 genes in the N. lugens genome. P450 sequences were obtained from the N. lugens genomic sequences and transcriptome databases, which were confirmed using the tBLASTX algorithm with a cut-off E-value of 10−10. The genomic organization of exons and introns of the P450s was predicted based on mRNA–genome alignments at the NCBI Spideyweb (http://www.ncbi.nlm.nih.gov/spidey/spideyweb.cgi). Locus, size and orientation indicate the location on scaffolds, predicted amino acids (aa) and the transcription orientation of the genes. L. striatellus, Laodelphax striatellus; S. furcifera, Sogatella furcifera; R. flavipes, Reticulitermes flavipes; L. m. manilensis, Locusta migratoria manilensis; A. pisum, Acyrthosiphon pisum; A. aegypti, Aedes aegypti; D. citri, Diaphorina citri; M. sexta, Manduca sexta; Z. nevadensis, Zootermopsis nevadensis. Predicted gene

GenBank ID

Locus

Size (aa)

Exon

Orientation

Best match

Similarity

CYP2 (10) NlCYP15G1 NlCYP15G1 NlCYP18A1 NlCYP303A1 NlCYP304H1v4 NlCYP304H1v4 NlCYP305A15 NlCYP306A2 NlCYP307A2 NlCYP307B1

KM217004 KM217006 KM217007 KM217008 KM217009 KM217010 KM217011 KM217013 KM217014 KM217015

scaffold1053 scaffold2132 scaffold953 scaffold1551 scaffold1903 scaffold448 scaffold1188 scaffold953 scaffold1880 scaffold838

495 495 532 503 506 506 494 499 520 507

7 2 7 11 8 13 9 6 6 3

− − + − − + − − − −

R. flavipes R. flavipes L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella S. furcifera L. striatella

65% 65% 98% 98% 95% 95% 91% 91% 91% 94%

CYP3 (19) NlCYP408A1v2 NlCYP6AX1v2 NlCYP6AX1v2 NlCYP6AY1v2 NlCYP6BD12 NlCYP6CS1v2 NlCYP6CW1v2 NlCYP6ER1v2 NlCYP6FK2 NlCYP6FL4 NlCYP6FL4 NlCYP6FL3 NlCYP6FJ2 NlCYP3115A1 NlCYP6FU1 NlCYP418A1v2 NlCYP427A1 NlCYP427A1 NlCYP427A1

KM217016 KM217017 KM217018 KM217019 KM217020 KM217021 KM217022 KM217023 KM217024 KM217025 KM217026 KM217027 KP710972 KM217028 KM217029 KM217031 KM217030 KM217032 KM217033

scaffold868 scaffold297 scaffold297 scaffold469 scaffold213 scaffold53 scaffold215 scaffold2255 scaffold53 scaffold2212 scaffold2990 scaffold700 scaffold1125 scaffold191 scaffold78 scaffold3261 scaffold3261 scaffold59 scaffold59

516 514 514 507 520 521 540 506 514 532 532 530 527 532 509 500 503 503 503

16 1 4 6 9 6 1 7 4 4 4 4 6 6 6 4 8 8 7

− + + − + + − + + + − − + + − − − + +

L. L. L. L. L. L. L. L. L. L. L. L. L. L. L. L. L. L. L.

59% 95% 95% 94% 94% 93% 91% 83% 92% 61% 61% 60% 90% 54% 93% 87% 93% 92% 92%

CYP4 (27) NlCYP4C77 NlCYP4C61v2 NlCYP4C76 NlCYP380C10 NlCYP4C78 NlCYP4C62v2 NlCYP4CE1v2 NlCYP4DE1 NlCYP4DE1 NlCYP4DC1 NlCYP4DD1 NlCYP4G76 NlCYP4G115 NlCYP4FB1 NlCYP4FB2 NlCYP417A1v2 NlCYP417A2 NlCYP417A3 NlCYP417B1 NlCYP417B1 NlCYP425A1 NlCYP425A1 NlCYP425B1 NlCYP426A1 NlCYP439A1 NlCYP439A2 NlCYP439B1

KM217034 KM217037 KM217038 KM217035 KM217036 KM217039 KM217040 KM217041 KM217042 KM217043 KM217044 KM217045 KM217046 KM217047 KM217048 KM217051 KM217050 KM217049 KM217052 KM217053 KM217055 KM217056 KM217057 KM217058 KM217059 KM217061 KM217060

scaffold745 scaffold745 scaffold745 scaffold559 scaffold1709 scaffold2554 scaffold232 scaffold132 scaffold960 scaffold922 scaffold107 scaffold3792 scaffold6648 scaffold439 scaffold895 scaffold371 scaffold371 scaffold371 scaffold926 scaffold926 scaffold2881 scaffold1253 scaffold1603 scaffold590 scaffold196 scaffold1840 scaffold10

509 505 512 547 521 520 541 534 534 502 523 572 578 293 509 510 507 505 515 515 513 513 522 255 400 354 437

12 12 12 16 13 11 14 4 9 9 7 5 6 4 7 11 11 10 12 12 5 2 3 3 5 2 9

+ + + + − − − − + + − + + + − − − + + + + + − − − − −

L. striatella L. striatella L. striatella A. pisum A. aegypti L. striatella Z. nevadensis L. striatella L. striatella L. striatella L. striatella D. citri M. sexta L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella L. striatella

96% 93% 90% 74% 79% 98% 68% 98% 98% 93% 97% 78% 78% 65% 61% 93% 99% 81% 97% 97% 97% 97% 94% 96% 96% 81% 67%

Mitochondrial P450 (12) NlCYP301A1 NlCYP301B1 NlCYP302A1 NlCYP302A1 NlCYP314A1 NlCYP314A1 NlCYP315A1 NlCYP353D1

KM216992 KM216993 KM216994 KM216995 KM216996 KM216997 KM216998 KM216999

scaffold287 scaffold644 scaffold63 scaffold3233 scaffold639 scaffold249 scaffold203 scaffold1093

521 532 563 563 542 542 480 472

3 10 6 6 7 3 8 4

− − − − − + − +

L. L. L. L. L. L. L. L.

98% 97% 94% 94% 97% 97% 93% 93%

m. manilensis striatella striatella striatella striatella striatella striatella striatella striatella striatella striatella striatella striatella striatella striatella striatella striatella striatella striatella

striatella striatella striatella striatella striatella striatella striatella striatella

(continued on next page)

Please cite this article as: S.-H. Lao, et al., Genomics (2015), http://dx.doi.org/10.1016/j.ygeno.2015.07.010

4

S.-H. Lao et al. / Genomics xxx (2015) xxx–xxx

Table 1 (continued) Predicted gene

GenBank ID

Locus

Size (aa)

Mitochondrial P450 (12) NlCYP353D1 NlCYP404B2 NlCYP404B2 NlCYP419A1

KM217000 KM217001 KM217002 KM217003

scaffold950 scaffold959 scaffold464 scaffold1908

472 481 481 658

samples. Eleven P450 genes exhibited notably increased expression levels upon imidacloprid stimulation (Fig. 4). These genes include a NlCYP18A1 in the CYP2 clade; NlCYP6FL4, NlCYP6ER1v2 and NlCYP6CS1v2 in the CYP3 clade; and NlCYP4DE1, NlCYP380C10, NlCYP4CE1v2, NlCYP4C76, NlCYP4FB1, NlCYP417A1v2 and NlCYP439A1 in the CYP4 clade. The expression level of NlCYP4CE1v2 increased more than sixfold and the other P450 gene expression levels increased approximately 2–3-fold after imidacloprid treatment. In contrast, only a few P450 genes showed significant changes in expression levels upon triazophos or deltamethrin treatment. The expression levels of the NlCYP303A1, NlCYP426A1 and NlCYP304H1v4 genes were up-regulated more than threefold by triazophos and/or deltamethrin, but were not affected by imidacloprid. Among these genes, NlCYP303A1 expression levels changed in response to triazophos and deltamethrin (Fig. 4). Triazophos or deltamethrin treatments led to significantly reduced expression levels of NlCYP6CS1v2. Intriguingly, the expression level of NlCYP6CS1v2 was up-regulated by imidacloprid (Fig. 4). The expression profile changes of N. lugens P450 genes in response to the different insecticides reflected that this insect species employs multiple defense strategies against a wide variety of environmental toxicants. We next examined the P450 gene expression patterns in insecticidetreated N. lugens compared to the acetone- and non-treated samples using quantitative real-time PCR. We confirmed that the expressions of five genes were significantly responsive to imidacloprid stimulation (Fig. 5). NlCYP6CS1v2 and NlCYP4CE1v2 expressions were strongly induced by imidacloprid, resulting in greater than threefold changes in expression levels compared to acetone and no treatment. Three

Exon 2 1 5

Orientation

Best match

Similarity

+ + + −

L. L. L. L.

93% 98% 98% 97%

striatella striatella striatella striatella

P450 genes (NlCYP4DE1, NlCYP417A1v2 and NlCYP439A1) had elevated expression levels, each notably increased compared to no treatment and more than twofold compared to acetone treatment. We also investigated the expression change of NlCYP6AY1v2 as it was previously implicated in imidacloprid resistance in this insect. RNA-Seq data showed increased expression (log2 ratio = 0.55) of NlCYP6AY1v2 in the 5th instar nymphs after imidacloprid treatment. Quantitative realtime PCR analysis confirmed that the expression levels increased 1.5fold compared to acetone and no treatment. NlCYP303A1 gene expression was dramatically up-regulated upon triazophos and deltamethrin treatments — transcript levels increased more than eightfold compared to no treatment and at least twofold compared to acetone treatment. 2.3. Development- and sex-specific expression analysis To understand the potentially physiological functions of N. lugens P450 genes, we analyzed the expression specificity in the different developmental stages and sexes using quantitative real-time PCR. These genes showed distinct expression patterns (Fig. 6). Transcripts of two mitochondrial P450 genes, (NlCYP301A1 and NlCYP301B1) were detected at the highest levels in eggs followed by 2nd and 5th instar nymphs, but were at very low levels in male and female adults. NlCYP306A2, NlCYP307B1 and NlCYP4G76 also displayed similar expression patterns. Their transcript levels gradually decreased between eggs and adults. RNA-Seq data showed that imidacloprid, triazophos and deltamethrin did not significantly change their expression levels, suggesting that these P450 genes did not respond to the insecticides; instead they

Fig. 1. The distribution of N. lugens P450 genes on scaffolds. The exon–intron organization of P450s was determined by matching the open reading frames of the transcript sequences to N. lugens genomic sequences with the online tool Spidey (http://www.ncbi.nlm.nih.gov/spidey/spideyweb.cgi). The orange arrows indicate the transcriptional orientation and sizes of P450s on scaffolds. Exons are indicated by green boxes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: S.-H. Lao, et al., Genomics (2015), http://dx.doi.org/10.1016/j.ygeno.2015.07.010

S.-H. Lao et al. / Genomics xxx (2015) xxx–xxx

5

Fig. 2. Phylogenetic analysis of insect P450 genes. The phylogenetic tree was constructed based on the predicted amino acid sequences of the P450 genes using Mega 5.05 (http://www. megasoftware.net/) with the Maximum Likelihood method. The Jones–Taylor–Thornton amino acid substitution model was used, while a test of phylogeny was carried out using a bootstrap analysis of 1000 replications; bootstrap values N50% are shown on each node of the tree. Green, pink, blue and red boxes represent the CYP2–CYP4 and mitochondrial P450 clades, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

might have physiological functions in embryonic development. Two members of the CYP2 clade (NlCYP303A1 and NlCYP304H1v4) showed a distinct expression pattern with notably high transcript levels in 2nd instar nymphs but low levels in 5th instar nymphs, adults and eggs, implying that they likely play important roles in early developmental stages. NlCYP303A1 transcript levels in 5th instar nymphs were strongly induced by triazophos and deltamethrin, but not by imidacloprid (Fig. 5), indicating potential sensitivity to organophosphates and pyrethroids. Interestingly, 10 P450 genes (NlCYP6FU1, NlCYP6FL4, NlCYP427A1, NlCYP4C77, NlCYP4CE1v2, NlCYP4DE1, NlCYP4G115, NlCYP439A1, NlCYP380C10 and NlCYP417B1) had common high expression in male adults. Among these, NlCYP4CE1v2, NlCYP4DE1, NlCYP439A1 and NlCYP417B1 exhibited much higher transcript levels in male compared to female adults. Their transcripts were detected at very low levels in eggs and nymphs. In contrast, NlCYP6AY1v2, NlCYP314A1 and NlCYP4C62v2 had significantly high expression in female adults. These P450s may be involved in the male- or femalerelated physiological processes. Among them, NlCYP4CE1v2, NlCYP4DE1, NlCYP439A1 and NlCYP6AY1v2 transcript levels in 5th instar nymphs were elevated by imidacloprid, but unaffected by triazophos and deltamethrin (Fig. 5), implying that they were sensitive or had potential resistance to imidacloprid treatment. The widely different expression patterns suggest that N. lugens P450s have multiple functions in biological processes in this insect species. The detailed functions and relationship between these P450 genes need to be further clarified.

3. Conclusions Completion of whole-genome sequencing of N. lugens enables the full identification of large and complex gene families, such as the P450 detoxification enzymes. Compared with studies that reported 32 fulllength or partial P450-related gene sequences from large-scale expressed sequence tag (EST) sequencing projects [26,27], the genome- and transcriptome-wide analysis in the present study identified a repertoire of P450 genes in N. lugens, which greatly extends our knowledge of the detoxification enzymes produced by this species. The available wholegenome sequences of insect species have shown that the most abundant P450 genes were in Culex quinquefasciatus (202), followed by A. aegypti (179) and Tribolium castaneum (131) [19]. In contrast, N. lugens had comparatively fewer P450 genes (68). Interestingly, the plant sap-sucking A. pisum encodes more P450 genes (83), consistent with its polyphagous habits and its exposure to a wider variety of plant secondary metabolites than N. lugens, thereby necessitating a larger complement of detoxification enzymes. The contraction of N. lugens P450 genes may be a consequence of its more specialized feeding habits, as it depends exclusively on rice phloem sap and would therefore require fewer detoxification enzymes than A. pisum. In N. lugens, the presence of CYP4CE1v2 and the absence of CYP6a2, CYP6a13, CYP6a14, CYP6k1 and CYP3A18 seemed to be specific to this insect species. A CYP6a2 gene has been shown to promote the metabolism of various insecticides in D. melanogaster [28]. The absence of CYP6a genes and a CYP3A18 gene in N. lugens might represent

Please cite this article as: S.-H. Lao, et al., Genomics (2015), http://dx.doi.org/10.1016/j.ygeno.2015.07.010

6

S.-H. Lao et al. / Genomics xxx (2015) xxx–xxx

Fig. 3. Gene Ontology (GO) classification of the differentially expressed genes upon insecticide treatments. The N. lugens genes responsive to insecticide treatments were annotated based on GO classification into three main categories: biological process, cellular component and molecular function. The Y-axes indicate the number of genes in a category. The upper panel (A) indicates the significantly differentially expressed genes between the imidacloprid- and acetone-treated samples; the lower panel (B) indicates the significantly differentially expressed genes between triazophos- and acetone-treated, and between deltamethrin- and acetone-treated samples. Control refers to acetone-treated samples.

an evolutionary strategy. Their lost functions have been most likely taken over by other P450 genes. The gene expression profiling of responses to insecticides provided new insights for understanding the potential functions of N. lugens P450s in physiological processes. Some P450 gene expressions were significantly increased by imidacloprid, but were reduced or unaffected by triazophos or deltamethrin; while several were notably enhanced by triazophos and/or deltamethrin, but not affected by imidacloprid, suggesting the sensitivity of the gene products to the different insecticides. The development- and sex-expression specificity provided useful information for a better understanding of the possible roles of N. lugens P450s in the reproduction and development processes. For example, in insects, Halloween genes that encode steroidogenic CYPs, including spook-CYP307A1, phantom-CYP306A1, disembodied-CYP302A1, shadow-CYP315A1 and shade-CYP314A1 regulated the pathway leading to the biosynthesis of the steroid hormone, 20hydroxyecdysone (20E). A step (the Black Box) in this pathway has not been fully characterized yet. It was thought that more than one steroidogenic CYP may participate in this step. Our results revealed that the homologs (NlCYP307B1 and NlCYP306A2) of spook-like and phantom-like genes had significantly high expressions in eggs, and a homolog (NlCYP314A1) of shade gene was highly expressed in female adults, suggesting that these P450s were most likely involved in embryonic development or reproduction. Several other genes including NlCYP301A1, NlCYP301B1, NlCYP4G76, NlCYP6AY1v2 and NlCYP4C62v2 displayed similar expression patterns to these Halloween-like genes. These findings provide valuable clues to explore the possible links to the steroid biosynthesis mechanisms during development processes. It is of interest and necessary to determine the detailed functions of P450 gene products in this insect species, which could produce potential targets for future management of this pest.

4. Materials and methods 4.1. Insects and insecticides The N. lugens populations used for RNA sequencing in this work were originally collected from a rice field in the Huajiachi Campus of Zhejiang University, Hangzhou, China, in 2008. One male and one female were mated to produce an F1 progeny. A single pair was then selected for sibling inbreeding for 13 generations. The purified colony was used for genomic DNA sequencing [19] and the experiments performed in this study. The insects were reared at 26 ± 0.5 °C with 50 ± 5% humidity on rice seedlings (Oryza sativa strain TN1) under a 16:8 hour light:dark photoperiod as previously described [20–22]. Technical imidacloprid, triazophos and deltamethrin were supplied by Xinnong Chemical Industrial Group Co. Ltd. (Taizhou, Zhejiang, China). 4.2. Generation of a digital gene expression library using high-throughput RNA-Seq The 5th instar N. lugens nymphs were used for RNA sequencing. Bioassays showed the 50% lethal doses (LD50) of triazophos, imidacloprid and deltamethrin (10.35, 7.72 and 7.38 ng/nymph, respectively). Under carbon dioxide anesthesia, each nymph was dosed with 0.25 μL of the insecticide solution on the prothorax notum using a hand microsyringe (ITO Standard Syringe, Shizuoka, Japan) mounted on a PB600-1 repeating dispenser (Hamilton Co, Reno, NV, USA). Acetonetreated nymphs were used as a control. After each insecticide treatment, nymphs were reared on fresh TN1 rice seedlings at 26 ± 0.5 °C with 50 ± 5% humidity under a 16:8 h light:dark photoperiod. Total RNA was extracted from the whole bodies of insects using Trizol reagent

Please cite this article as: S.-H. Lao, et al., Genomics (2015), http://dx.doi.org/10.1016/j.ygeno.2015.07.010

S.-H. Lao et al. / Genomics xxx (2015) xxx–xxx

7

Fig. 4. The N. lugens P450 gene expression changes in response to insecticide treatments. The 5th instar nymphs were treated with sub-lethal doses of imidacloprid, triazophos or deltamethrin. Total RNA was extracted from the whole body of each insecticide-treated N. lugens insect to construct the RNA-Seq cDNA libraries. Acetone-treated samples were used as controls. The differentially expressed P450s were identified based on the criteria of FDR ≤ 0.001 and log2 ratio ≥ 1 between each insecticide- and acetone-treated sample. Red, orange and yellow indicate genes with significantly increased expression levels, and blue and black indicate genes with significantly reduced expression levels. The Venn diagram (http://bioinfogp.cnb.csic. es/tools/venny/index.html) indicates the number of genes that were significantly differentially expressed in response to insecticide treatments. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

(Invitrogen, Carlsbad, CA, USA) at 48 h after treatment. In parallel, four digital gene expression libraries from acetone-, imidacloprid-, triazophos- and deltamethrin-treated samples were constructed using RNA-Seq. Briefly, mRNA was enriched from each total RNA sample using oligo (dT) magnetic beads and was interrupted to short fragments of approximately 200 nucleotides (nt). The first-strand cDNA was synthesized using random hexamer primers with the mRNA fragments as templates. After the second-strand cDNA synthesis, double-strand cDNA was purified with a QiaQuick PCR extraction kit and washed with EB buffer for end-repair and single adenine nucleotide addition. Finally, sequencing adaptors were ligated to the cDNA fragments. The cDNA was size-fractionated by agarose gel electrophoresis and enriched by PCR amplification. The library products were sequenced using the Illumina HiSeq™ 2000 platform. 4.3. Identification of P450 genes from the N. lugens genomic and transcriptomic databases P450 sequences were searched against the N. lugens genomic and transcriptomic databases based on the gene annotations as previously described [29]. To avoid possibly incomplete annotations, P450 sequences from other insect species [19] were used as queries in BLASTX searches against the N. lugens databases. The putative P450

genes were validated by aligning to the non-redundant (nr) National Center for Biotechnology Information (NCBI) nucleotide database using a cut-off E-value of 10− 10. The exon–intron organizations of P450 genes were predicted based on the mRNA–genome alignments at the NCBI Spideyweb (http://www.ncbi.nlm.nih.gov/spidey/ spideyweb.cgi).

4.4. Analysis of the expressed P450 genes responsive to insecticide treatments The N. lugens genomic and transcriptomic databases were used as reference sequences to map and analyze the insecticide-induced gene expression profiling. Prior to mapping and analysis, all sequencing reads were filtered to remove those reads that had only adaptors, N10% unknown (N) bases or low quality reads (i.e. the percentage of low quality bases with a quality value ≤ 5 was N50% in a read). The clean reads were mapped to the N. lugens reference sequence databases. The number of clean reads was calculated and then normalized to the TPM (number of transcripts per million tags). The differentially expressed P450 genes were identified based on the following criteria: FDR ≤ 0.001 and log2 ratio ≥ 1 between insecticide- and acetonetreated samples.

Please cite this article as: S.-H. Lao, et al., Genomics (2015), http://dx.doi.org/10.1016/j.ygeno.2015.07.010

8

S.-H. Lao et al. / Genomics xxx (2015) xxx–xxx

Fig. 5. Analysis of P450 gene expression in response to insecticide treatment. The 5th instar nymphs were treated with a sub-lethal dose of imidacloprid, triazophos or deltamethrin. Total RNA was extracted from the nymphs (n = 15–20) at 48 h after treatments. No- and acetone-treatment samples were used as controls. First-strand cDNA was analyzed in each quantitative real-time PCR reaction. The reactions were performed with specific primers to amplify the P450s. The relative expression levels of each gene were normalized using the N. lugens 18s rRNA threshold cycle (Ct) values. Three biological replicates were carried out for each gene based on independent samples and the relative transcript levels were calculated using the ΔΔCt method. Results of triplicate experiments are shown with the standard deviations.

4.5. Phylogenetic analysis The deduced N. lugens P450 amino acid sequences were aligned with the best-matched homologs of other insect species using the ClustalX

program [30]. The phylogenetic tree was generated using the Mega 5.05 program (http://www.megasoftware.net/) by the maximum likelihood method. The evolutionary relationships were determined by bootstrap analysis with 1000 replications.

Fig. 6. Analysis of development- and sex-expression specificity of N. lugens P450 genes. Total RNAs were extracted from different developmental stages and sexes of N. lugens (n = 10–80) and used for investigating P450 gene expression using quantitative real-time PCR analysis, as described in Fig. 5. The reactions were performed with specific primers for amplifying P450 genes. Three biological replications were carried out based on independent RNA sample. 2nd, 5th, F and M refer to 2nd and 5th instar nymphs and female and male adults, respectively. Red, green, blue and pink boxes indicate the overall common patterns of P450 genes that showed high expression in eggs, 2nd instar nymphs and female and male adults, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: S.-H. Lao, et al., Genomics (2015), http://dx.doi.org/10.1016/j.ygeno.2015.07.010

S.-H. Lao et al. / Genomics xxx (2015) xxx–xxx

4.6. Confirmation of differential expressions using quantitative real-time PCR Total RNA was extracted from N. lugens specimens using Trizol reagent (Invitrogen). RNA was reverse-transcribed using the ReverTra Ace® qPCR RT Master Mix with a gDNA Remover Kit (ToYoBo, Osaka, Japan) to remove any contaminating genomic DNA. Quantitative realtime PCR was run on a CFX96™ Real-Time System (Bio-Rad, Hercules, CA, USA) using the iQ™ SYBR Green® Supermix Kit (Bio-Rad). The first-strand cDNA and a no-reverse-transcription control were used as templates for three biological replicates under the following reaction program: an initial denaturation step at 95 °C for 3 min followed by 40 cycles at 95 °C for 5 s and 60 °C for 35 s. The gene-specific primers were designed using the Primer Premier 5 program based on the N. lugens transcriptomic sequences as shown in supplementary file Table S1. The N. lugens housekeeping gene for 18S ribosomal RNA was used as an internal control. The ΔΔCt method was used to evaluate the quantitative variation in the transcript levels as described previously [22,29,31]. 4.7. Availability of supporting data The supporting data in this study have been submitted to the open access repositories. N. lugens genome assemblies have been deposited at GenBank under accession number AOSB00000000 (BioProject PRJNA177647). The N. lugens transcriptomic dataset is available in the Sequence Read Archive (SRA) database (http://www.ncbi.nlm.nih.gov/ sra). The accession number of the N. lugens transcriptomic dataset is SRX023419. The transcript sequences of N. lugens P450 genes were submitted to the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/). The accession numbers of these sequences are as follows: KM216992–KM217061 and KP710972. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.ygeno.2015.07.010. Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant No. 31371934). References [1] Z. Wei, W. Hu, Q. Lin, X. Cheng, M. Tong, L. Zhu, R. Chen, G. He, Understanding rice plant resistance to the Brown Planthopper (Nilaparvata lugens): a proteomic approach, Proteomics 9 (2009) 2798–2808. [2] K.K. Jena, S.M. Kim, Current status of brown planthopper (BPH) resistance and genetics, Rice 3 (2010) 161–171. [3] C. Bass, R. Carvalho, L. Oliphant, A. Puinean, L. Field, R. Nauen, M. Williamson, G. Moores, K. Gorman, Overexpression of a cytochrome P450 monooxygenase, CYP6ER1, is associated with resistance to imidacloprid in the brown planthopper, Nilaparvata lugens, Insect Mol. Biol. 20 (2011) 763–773. [4] E. Cohen, Pesticide-mediated homeostatic modulation in arthropods, Pestic. Biochem. Physiol. 85 (2006) 21–27. [5] Y.H. Wang, S.G. Wu, Y.C. Zhu, J. Chen, F.Y. Liu, X.P. Zhao, Q. Wang, Z. Li, X.P. Bo, J.L. Shen, Dynamics of imidacloprid resistance and cross‐resistance in the brown planthopper, Nilaparvata lugens, Entomol. Exp. Appl. 131 (2009) 20–29. [6] K. Gorman, Z. Liu, I. Denholm, K.U. Brüggen, R. Nauen, Neonicotinoid resistance in rice brown planthopper, Nilaparvata lugens, Pest Manag. Sci. 64 (2008) 1122–1125. [7] Y. Wen, Z. Liu, H. Bao, Z. Han, Imidacloprid resistance and its mechanisms in field populations of brown planthopper, Nilaparvata lugens Stål in China, Pestic. Biochem. Physiol. 94 (2009) 36–42.

9

[8] L. Zewen, H. Zhaojun, W. Yinchang, Z. Lingchun, Z. Hongwei, L. Chengjun, Selection for imidacloprid resistance in Nilaparvata lugens: cross‐resistance patterns and possible mechanisms, Pest Manag. Sci. 59 (2003) 1355–1359. [9] A. Puinean, I. Denholm, N. Millar, R. Nauen, M. Williamson, Characterisation of imidacloprid resistance mechanisms in the brown planthopper, Nilaparvata lugens Stål (Hemiptera: Delphacidae), Pestic. Biochem. Physiol. 97 (2010) 129–132. [10] Z. Ding, Y. Wen, B. Yang, Y. Zhang, S. Liu, Z. Liu, Z. Han, Biochemical mechanisms of imidacloprid resistance in Nilaparvata lugens: over-expression of cytochrome P450 CYP6AY1, Insect Biochem. Mol. Biol. 43 (2013) 1021–1027. [11] X.C. Li, M.A. Schuler, M.R. Berenbaum, Molecular mechanisms of metabolic resistance to synthetic and natural xenobiotics, Annu. Rev. Entomol. 52 (2007) 231–253. [12] R. Feyereisen, Evolution of insect P450, Biochem. Soc. Trans. 34 (2006) 1252–1255. [13] M. Maïbèche-Coisne, A.A. Nikonov, Y. Ishida, E. Jacquin-Joly, W.S. Leal, Pheromone anosmia in a scarab beetle induced by in vivo inhibition of a pheromonedegrading enzyme, Proc. Natl. Acad. Sci. 101 (2004) 11459–11464. [14] R. Feyereisen, Insect P450 enzymes, Annu. Rev. Entomol. 44 (1999) 507–533. [15] K.F. Rewitz, L.I. Gilbert, Daphnia Halloween genes that encode cytochrome P450s mediating the synthesis of the arthropod molting hormone: evolutionary implications, BMC Evol. Biol. 8 (2008) 60. [16] R. Niwa, T. Matsuda, T. Yoshiyama, T. Namiki, K. Mita, Y. Fujimoto, H. Kataoka, CYP306A1, a cytochrome P450 enzyme, is essential for ecdysteroid biosynthesis in the prothoracic glands of Bombyx and Drosophila, J. Biol. Chem. 279 (2004) 35942–35949. [17] W. Mao, M.A. Berhow, A.R. Zangerl, J. Mcgovern, M.R. Berenbaum, Cytochrome P450-mediated metabolism of xanthotoxin by Papilio multicaudatus, J. Chem. Ecol. 32 (2006) 523–536. [18] J.S. Ramsey, D.S. Rider, T.K. Walsh, M. De Vos, K. Gordon, L. Ponnala, S. Macmil, B. Roe, G. Jander, Comparative analysis of detoxification enzymes in Acyrthosiphon pisum and Myzus persicae, Insect Mol. Biol. 19 (2010) 155–164. [19] J. Xue, X. Zhou, C.X. Zhang, L.L. Yu, H.W. Fan, Z. Wang, H.J. Xu, Y. Xi, Z.R. Zhu, W.W. Zhou, P.L. Pan, B.L. Li, J.K. Colbourne, H. Noda, Y. Suetsugu, T. Kobayashi, Y. Zheng, S.L. Liu, R. Zhang, Y. Liu, Y.D. Luo, D.M. Fang, Y. Chen, D.L. Zhan, X.D. Lv, Y. Cai, Z.B. Wang, H.J. Huang, R.L. Cheng, X.C. Zhang, Y.H. Lou, B. Yu, J.C. Zhuo, Y.X. Ye, W.Q. Zhang, Z.C. Shen, H.M. Yang, J. Wang, J. Wang, Y.Y. Bao, J.A. Cheng, Genomes of the rice pest brown planthopper and its endosymbionts reveal complex complementary contributions for host adaptation, Genome Biol. 15 (2014) 521. [20] J. Xue, Y.Y. Bao, B.L. Li, Y.B. Cheng, Z.Y. Peng, H. Liu, H.J. Xu, Z.R. Zhu, Y.G. Lou, J.A. Cheng, Transcriptome analysis of the brown planthopper Nilaparvata lugens, PLoS One 5 (2010) e14233. [21] Y.Y. Bao, Y. Wang, W.J. Wu, D. Zhao, J. Xue, B.Q. Zhang, Z.C. Shen, C.X. Zhang, De novo intestine-specific transcriptome of the brown planthopper Nilaparvata lugens revealed potential functions in digestion, detoxification and immune response, Genomics 99 (2012) 256–264. [22] Y.Y. Bao, L.Y. Qu, D. Zhao, L.B. Chen, H.Y. Jin, L.M. Xu, J.A. Cheng, C.X. Zhang, The genome- and transcriptome-wide analysis of innate immunity in the brown planthopper, Nilaparvata lugens, BMC Genomics 14 (2013) 160. [23] T. Namiki, R. Niwa, T. Sakudoh, K. Shirai, H. Takeuchi, H. Kataoka, Cytochrome P450 CYP307A1/Spook: a regulator for ecdysone synthesis in insects, Biochem. Biophys. Res. Commun. 337 (2005) 367–374. [24] C. Helvig, J. Koener, G. Unnithan, R. Feyereisen, CYP15A1, the cytochrome P450 that catalyzes epoxidation of methyl farnesoate to juvenile hormone III in cockroach corpora allata, Proc. Natl. Acad. Sci. 101 (2004) 4024–4029. [25] The International Aphid Genomics Consortium, Genome sequence of the pea aphid Acyrthosiphon pisum, PLoS Biol. 8 (2010) e1000313. [26] H. Noda, S. Kawai, Y. Koizumi, K. Matsui, Q. Zhang, S. Furukawa, M. Shimomura, K. Mita, Annotated ESTs from various tissues of the brown planthopper Nilaparvata lugens: a genomic resource for studying agricultural pests, BMC Genomics 9 (2008) 117. [27] C. Bass, M.B. Hebsgaard, J. Hughes, Genomic resources for the brown planthopper, Nilaparvata lugens: transcriptome pyrosequencing and microarray design, Insect Sci. 19 (2012) 1–12. [28] G. Goff, S. Boundy, P. Daborn, J. Yen, L. Sofer, R. Lind, C. Sabourault, L. Madi-Ravazzi, Microarray analysis of cytochrome P450 mediated insecticide resistance in Drosophila, Insect Biochem. Mol. Biol. 33 (2003) 701–708. [29] Y.Y. Bao, X. Qin, B. Yu, L.B. Chen, Z.C. Wang, C.X. Zhang, Genomic insights into the serine protease gene family and expression profile analysis in the planthopper, Nilaparvata lugens, BMC Genomics 15 (2014) 507. [30] J.D. Thompson, T.J. Gibson, F. Plewniak, F. Jeanmougin, D.G. Higgins, The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools, Nucleic Acids Res. 25 (1997) 4876–4882. [31] K.J. Livak, T.D. Schmittgen, Analysis of relative gene expression data using real-time quantitative PCR and the 2(− Delta Delta C(T)) method, Methods 25 (2001) 402–408.

Please cite this article as: S.-H. Lao, et al., Genomics (2015), http://dx.doi.org/10.1016/j.ygeno.2015.07.010