RNA-Seq and molecular docking reveal multi-level pesticide resistance in the bed bug

Background Bed bugs (Cimex lectularius) are hematophagous nocturnal parasites of humans that have attained high impact status due to their worldwide resurgence. The sudden and rampant resurgence of C. lectularius has been attributed to numerous factors including frequent international travel, narrower pest management practices, and insecticide resistance. Results We performed a next-generation RNA sequencing (RNA-Seq) experiment to find differentially expressed genes between pesticide-resistant (PR) and pesticide-susceptible (PS) strains of C. lectularius. A reference transcriptome database of 51,492 expressed sequence tags (ESTs) was created by combining the databases derived from de novo assembled mRNA-Seq tags (30,404 ESTs) and our previous 454 pyrosequenced database (21,088 ESTs). The two-way GLMseq analysis revealed ~15,000 highly significant differentially expressed ESTs between the PR and PS strains. Among the top 5,000 differentially expressed ESTs, 109 putative defense genes (cuticular proteins, cytochrome P450s, antioxidant genes, ABC transporters, glutathione S-transferases, carboxylesterases and acetyl cholinesterase) involved in penetration resistance and metabolic resistance were identified. Tissue and development-specific expression of P450 CYP3 clan members showed high mRNA levels in the cuticle, Malpighian tubules, and midgut; and in early instar nymphs, respectively. Lastly, molecular modeling and docking of a candidate cytochrome P450 (CYP397A1V2) revealed the flexibility of the deduced protein to metabolize a broad range of insecticide substrates including DDT, deltamethrin, permethrin, and imidacloprid. Conclusions We developed significant molecular resources for C. lectularius putatively involved in metabolic resistance as well as those participating in other modes of insecticide resistance. RNA-Seq profiles of PR strains combined with tissue-specific profiles and molecular docking revealed multi-level insecticide resistance in C. lectularius. Future research that is targeted towards RNA interference (RNAi) on the identified metabolic targets such as cytochrome P450s and cuticular proteins could lay the foundation for a better understanding of the genetic basis of insecticide resistance in C. lectularius.

It has been more than half a century since C. lectularius first showed resistance to DDT [11,12]. The sudden resurgence of C. lectularius is purportedly due to increased resistance to broad-spectrum insecticides, changes in pest management practices, frequent international travel and passive dispersal (clothing, luggage and second-hand furniture) [13,14]. In the recent past, C. lectularius severely affected the hospitality industry, wherein some hotels were closed due to heavy infestations [15]. The lack of effective pest management tools for this blood-feeding insect has led to its successful establishment around the globe [8].
Insects develop resistance to insecticides through four modes: penetration resistance (thicker cuticle for decreased entry of insecticides), behavioral resistance (avoidance of the toxic compounds), target site resistance (knockdown resistance, kdr), and metabolic resistance (detoxification primarily through the action of cytochrome P450s, glutathione S-transferases [GSTs], and carboxylesterases) [16]. Among these four modes, target site resistance (kdr) has been well characterized in C. lectularius [17,18], but little to no knowledge exists regarding the other modes of resistance. Recent transcriptomic studies of C. lectularius have set the foundation for understanding the potential contribution of cytochrome P450s (metabolic resistance) in pesticide resistance [19,20]. However, a more comprehensive and/ or global understanding of the genes involved in pesticide resistance and their "modus operandi" in C. lectularius are necessary for understanding the genetic factors establishing resistance in C. lectularius and for improving existing control strategies or devising new ones. The next-generation sequencing (NGS) methods via different platforms (Illumina Genome Analyzer, Applied Biosystems SOLiD, Helicos Biosciences Heliscope, Roche 454 Life Sciences) have revolutionized functional genomics research in non-model organisms [21][22][23][24]. Illumina deep sequencing (RNA sequencing, RNA-Seq) has emerged as a powerful tool for simultaneous transcriptome characterization and differential gene expression (DGE) analysis to better understand eco-physiological adaptations of insects [24][25][26][27][28].
In this study, we used RNA-Seq to compare and contrast pesticide-resistant (PR) and pesticide-susceptible (PS) strains of C. lectularius. As a result, we have significantly enriched the existing transcriptomic database of C. lectularius [19] and subsequently identified candidate genes putatively involved in insecticide resistance. Lastly, molecular docking was performed for a candidate CYP3 clan member (CYP397A1V2) to further determine the likely contribution of P450 mediated insecticide resistance in C. lectularius.

Results and Discussion
De novo assembly We generated 62,107,336 and 64,214,910 reads for two PR strains and 72,748,924 and 72,400,340 reads for two PS strains (Table 1). An overview of the de novo assembly is shown in Figure 1. De novo assembly of the transcriptomes using an automated transcriptome assembly pipeline (see Methods) with 'PR' sample 1 and ''PS' sample 1 resulted in 34,385 (with N50 being 833 nt) and 46,412 (with N50 being 1064 nt) contigs, respectively (Table 1). Combining these contigs with previously assembled contigs (21,088; with N50 being 466 nt) [19] resulted in a total of 51,492 ESTs (after removal of redundant sequences) and these were used as reference sequences (Table 1, Figure 2A). The ESTcalc http://fgp. huck.psu.edu/NG_Sims/ngsim.pl estimated that the combined Illumina reads and 454 datasets would cover 90% of the bed bug transcriptome (Additional file 1).
The read variation number observed in the current study has been commonly associated with the Illumina platform as different lanes can produce different number of counts [25,27]. Other possible reasons for read variation include different percentages of poly A + in different samples, differential efficiency of conversion of RNA to cDNA, size or GC composition variation among libraries (smaller and/or lower GC content amplicons amplify more efficiently during PCR for library quantitation). All these factors may contribute to slight variations in the estimation of the library concentration that is applied onto the flow-cell, finally resulting in differences in cluster densities. The coverage and number of reads (i.e., depth of sequencing) obtained in this study far exceeds that of another recently published article on the C. lectularius transcriptome involving differential expression between susceptible and resistant populations [20].
To determine the completeness of our final assembly, each of the derived assembled sequences was compared with its putative ortholog from Acyrthosiphon pisum, Anopheles gambiae, Drosophila melanogaster, and Pediculus humanus corporis. Nearly, 33-35% of the assembled sequences exceeded a 0.8 ortholog hit ratio and 48-50% had more than a 0.5 ortholog hit ratio (Figure 2B). An ortholog ratio of one indicates a fully assembled transcriptome, while a value close to zero indicates a poor assembly [29]. The current transcriptome characterization resulted in a 23-fold enrichment of the existing EST database of C. lectularius (1,971 ESTs http://www.ncbi.nlm.nih.gov accessed on 25 th August 2011) which allowed fine characterization of the C. lectularius transcriptome as has been observed in many non-model insect species [24][25][26][27][28].

Comparative genomics
About 43% (21,908/51,492) of the C. lectularius ESTs had one or more hits to protein sequences in the nonredundant (nr) protein database http://www.ncbi.nlm. nih.gov with the remaining sequences (57%) being transcripts of unknown function (TUF). This is in agreement with our previous 454 EST datasets of C. lectularius, wherein we reported a similar percentage (45.2%) of EST hits with the nr protein database [19]. The majority of TUF sequences might be due to novel transcripts or genes whose biological functions are not assigned as has been observed in other transcriptomic studies [19,[25][26][27]. The top BLAST hits of the C. lectularius ESTs showed majority hits to insects (65.45%) and non-insect eukaryotes (32.82%); there were a few hits (1.73%) to members of bacteria (Additional file 2). This trend in top BLAST hits is consistent with other NGS datasets reported for insects [19,30,31]. Similarity searches between C. lectularius ESTs and protein sequences of Ac. pisum, A. gambiae, D. melanogaster and P. humanus corporis revealed an overall 41% similarity (21,157/51,492), with the highest number of hits to Ac. pisum (36%; 18,568/51,492) and P. humanus corporis (35.8%; 18,456/51,492) ( Figure 3). This similarity of C. lectularius with Ac. pisum and P. humanus corporis might be due to their phylogenetic relatedness.

Functional annotation
A total of 15,540 out of the 51,492 ESTs (reference sequences) were annotated and Gene Ontology (GO) terms were found to be distributed in a wide variety of functional categories (Additional file 3). Annotation with the D. melanogaster genome revealed a similar distribution  of GO categories with C. lectularius showing no notable bias towards any category. The majority of the ESTs annotated with biological processes were involved in development and biological regulation, while the ESTs annotated with molecular function revealed high catalytic and binding activity (Additional file 3). A number of putative viral sequences were represented in the GO terms for biological processes (viral life cycle) and cellular component (viroids) (Additional file 3 and Additional file 4). Though viral transmission by C. lectularius is still not clear, the Hepatitis B virus (HBV) has frequently been detected in wild C. lectularius [7]. Future studies are required to both validate these viral sequences and determine their biological relevance within the C. lectularius genome. The assigned metabolic pathways revealed a high number of sequences (1,557/8,159) to be involved in synthesis of carbohydrates, proteins, lipids and nucleotides (Additional file 5). These predicted pathways together with gene annotation will better help in revealing gene function in C. lectularius [19].

Differential gene expression
RNA-Seq libraries in the current study were constructed from paired biological replicates of resistant and susceptible RNA samples of C. lectularius. The four RNA samples were processed in two batches: Batch-1 (Resistant-1, Susceptible-1) and Batch-2 (Resistant-2, Susceptible -2). Hence, we applied GLMseq (a custom R-script for fitting a two-way Generalized Linear Model) analysis to the four un-normalized RNA-Seq counts at each EST. To find differentially expressed ESTs, Illumina counts of all four populations (PS1, PS2, PR1and PR2, refer to the methods section for more details) were aligned to the assembled reference database which resulted in~15,000 highly significant differentially expressed ESTs (P < 0.005). All differentially expressed ESTs, along with their Absolute Log2 Fold Change Adjusted (ALFCA) and their description based on BLAST hits, are presented in Additional file 6. The primary objective of GLMseq analysis was to determine the amount and significance of the main resistance effects, setting aside any secondary batch effects due to transcript-specific interactions between batch effects and resistance effects. Currently available methods for RNA-Seq data analysis, such as DESeq, are limited to one-way analysis of replicate groups and consequently cannot correctly handle two-way interactions. Therefore, a custom R-script (GLMseq) was used to fit a two-way generalized linear model consisting of a main resistance effect and a secondary interaction effect.

GO analysis of differentially expressed ESTs
GO enrichment analysis was performed for the top 5,000 differentially expressed ESTs obtained in GLMseq and all of the molecular function, biological process and cellular components are shown in Additional file 7. We focused on molecular function in order to identify potential categories of genes that are associated with pesticide resistance, such as have been revealed in other insect studies [32]. Of the top 5,000 genes about 52.62% (2,631) were up-regulated and 47.38% (2,369) were down-regulated in the PR strains. Among these, 19.86% (993/5,000) of the ESTs had no known function (transcripts of unknown function; TUFs) (Additional file 7). For the up-regulated cluster of differentially expressed ESTs, some of the enriched GO terms for molecular function included ATP binding, actin binding, serine/ threonine kinase activity, structural constituent of cuticle, etc. ( Figure 4A). The enriched GO terms for actin binding in PR strains of C. lectularius might suggest their potential involvement in cytoskeletal networks and cellular mechanical integrity as observed in other pyrethroid-resistant insect strains [33].

Detoxification enrichment
We found enriched ATP binding and glutathione peroxidase-associated GO terms in the up-regulated cluster of differentially expressed ESTs suggesting the involvement of other detoxification genes such as ATP binding cassette (ABC) transporters and quenchers of reactive oxygen species (ROS) ( Figure 4A), which is in agreement with other insect studies [34,35]. High occurrences of P450s (29) in the top 5,000 differentially expressed ESTs of the C. lectularius PR strain suggest their putative role in metabolic resistance ( Table 2). Insect P450s are one of the key players in detoxification and metabolism of a broad range of toxins including plant-synthesized and synthetic compounds [36]. The top differentially expressed ESTs were found in all the major clades of P450s (25 CYP3, 2 CYP4, and 2 mitochondrial CYP clan members) (Additional file 8). Among these, the up-regulated differentials in PR strains (65.51%) were exclusively CYP3 clan members, which are thought to be the primary P450s involved in detoxification of xeno-and endobiotics in insects [37]. Besides the CYP3 clan of P450s, CYP4 and mitochondrial CYP clan members are also thought to participate in insecticide resistance [38]. However, in the top 5,000 differentials, the CYP4 and mitochondrial P450 members were found to be downregulated (Additional file 8).

Cuticular enrichment
The enriched GO terms for actin binding, structural constitution of cuticle, and structural constituent of muscle in PR strains of C. lectularius could imply the strengthening of structural components such as the cuticle and midgut-associated structures. Indeed among the up-regulated cluster of ESTs, we found a high number (46) of transcripts encoding for cuticular proteins ( Table 2). Cuticular proteins are major components of insect cuticle. The constitutive or induced expression of such proteins potentially establishes penetration resistance in C. lectularius, as observed in other insect systems [39]. It is interesting to observe that only 8 transcripts encoding cuticular proteins were found in the down-regulated differential cluster in PR strains of C. lectularius (Table 2).
For the down-regulated cluster of differentially expressed ESTs, the most enriched GO terms included metalloexopeptidase, aminopeptidase activities, manganese ion binding, lipid transporter activity, etc. ( Figure  4B). These trends could imply the transcriptomic adjustments within C. lectularius upon encountering insecticides. While Adelman et al. [20] attribute multilevel insecticide resistance of bed bugs primarily toward the activity of metabolic genes (cytochrome P450, carboxylesterases and GSTs) plus kdr mutations, their analysis failed to identify major players involved in penetration resistance (e.g., cuticular proteins). Further, in addition to the metabolic genes identified in their study, our unbiased/global approach identified other key players of metabolic resistance such as ABC transporters and antioxidant genes ( Table 2). We also screened the ESTs that were uniquely expressed in susceptible and resistant populations of C. lectularius, wherein a high number of ESTs (77) were found in PS strains and a few (7 ESTs) in PR strains (Additional File 9). In both cases, a majority of the ESTs were TUFs (85.71% in PR and 74.02% in PS strains). At the current time, it is too speculative to explain the uniquely expressed ESTs encoding for cathepsin in PR strains and RP45 and keratin-associated proteins among PS strains (Additional file 9).

RNA-Seq validation of candidate genes
To validate the expression profiles obtained through GLMseq analysis, quantitative real time RT-PCR (qRT-PCR) was performed on 12 selected candidate genes (cuticular proteins, CYP3 and mitochondrial P450 clan members, ABC transporter, superoxide dismutase, and acetylcholinesterase) belonging to the top differential cluster (Additional file 10). Except for one of the candidate genes profiled (mitochondrial CYP member, CYP 301A2), the qRT-PCR results correlated with the GLMseq profiles (Additional file 10). The expression patterns of CYP301A2 of C. lectularius are in corroboration with DDT-resistant and DDT-susceptible fruit flies (D. melanogaster), wherein no appreciable differences of CYP301A1 transcript levels were reported between the susceptible and resistant flies. Further, DDT treatment did not induce the expression levels of CYP301A1 [40]. The mitochondrial P450 clan members are exclusively found in animals with conserved and diversified groups potentially involved in physiological and detoxification processes in insects [41,42]. Though these clan members have been well documented to detoxify insecticides and plant allelochemicals [42][43][44], the CYP301A2 of C. lectularius might plausibly also participate in the synthesis of 20-hydroxyecdysone (20-HE) as observed in other insect studies [45,46]. ) and are consistent with the RNA-Seq profiles (Additional file 10). The higher expression of LCP and PCP in PR strains could contribute to increased pyrethroid resistance [47]. In addition to their role in other physiological processes, insects commonly up-regulate transcript levels of cuticular proteins to reduce the penetration rate of insecticides [39,47,48]. Higher transcript levels for CDA and CHS could contribute altered chitin and chitosan ratios in cuticle, which influences insect survival [49]. Co-expression of the CHS and CDA transcripts might potentially reduce the entry of xenobiotic compounds into the insect body. The higher expression of CPAP in PR strains (Additional file 10) indicates a possible role of the peritrophic membrane (gut) in sequestering, detoxifying ingested xenobiotics in addition to its role in peritrophic membrane formation and protection from invasive parasites [50,51].
All three cytochrome P450s of CYP3 clan members (CYP397A1V2, putative CYP6A2 and CYP6A13) revealed higher mRNA levels in the PR strains (Additional file 10). This is in agreement with other insect systems showing resistance to DDT or pyrethroids [52]. The transcript levels of a superoxide dismutase and an ABC transporter were also found to be higher in the PR strains (Additional file 10), which suggests that these proteins are potentially involved in the elimination and efflux of intracellular toxins thereby reducing their interaction with intracellular targets [35,53]. The transcripts encoding acetylcholinesterase (insecticide target protein) showed significantly higher expression in the PR strains (Additional file 10) possibly resulting in increased production of acetylcholinesterase to minimize the toxic effects of insecticides [54]. The current qRT-PCR profiles of three P450s and an ABC transporter together with our previous findings of a Phase II protein (GST) [19] suggest a phase-wise detoxification of xenobiotic compounds in C. lectularius. Further, the qRT-PCR profiles also validate our GLMseq analysis employed for generating the differentially expressed ESTs.
Tissue and development-specific expression of CYP397A1V2 and CYP6A13 Given our interest in the role of P450s in metabolic resistance, we further profiled the transcript levels of CYP397A1V2 and CYP6A13 in different tissues (cuticle, midgut and Malpighian tubules) of adults and different developmental stages (early and late instar nymphs and adults) of the PR and PS strains. Both genes showed lesser transcript levels in PS strains, therefore the PS strains were used as a calibrator to determine fold change in tissues and developmental stages of the PR samples.
An interesting finding of this study was the peak mRNA levels of CYP397A1V2 in the cuticle ( Figure 5A) and relatively high mRNA levels in the Malpighian tubules and midgut. This suggests the possible function of P450-mediated detoxification of pesticides in the cuticle of C. lectularius and/or its participation in other physiological functions. However, these findings need further functional validation. High transcript levels of both the P450s (CYP397A1V2 and CYP6A13) in the Malpighian tubules could also correlate with detoxification given their role in metabolism and excretion of endogenous solutes and xenobiotics [55]. The role of P450s in the insect midgut has been well established in detoxification as well as possible pheromone synthesis [56].
The qRT-PCR analysis of CYP397A1V2 and CYP6A13 showed higher mRNA levels in all developmental stages of PR compared to PS strains, specifically peak levels in early instar nymphs for CYP397A1V2 and late instar nymphs for CYP6A13 ( Figure 5B). Our findings correlate with those in other insect systems and therefore suggest the up-regulation of P450s in C. lectularius to efficiently detoxify a broad range of toxic substrates at different tissue interfaces and during development [19,42,[56][57][58].
The expression of CYP397A1V2 transcripts in adults are in agreement with recent transcriptomic studies of C. lectularius [20]. However, our results clearly demonstrate the broader appearance of key metabolic players of both penetration resistance as well as metabolic resistance compared to Adelman et al. [20].

Homology modeling
Given the unique tissue-specific expression profiles of CYP397A1V2 in PR strains, the full-length sequence of CYP397A1V2 was obtained, which revealed all the characteristic features of a CYP3 clan member (Additional file 11). A three-dimensional model for CYP397A1V2 was generated using human cytochrome P450 CYP3A4 (PDB code: 1WOE) as the template (Figure 6A and Additional file 12). During the homology modeling process, an iron-heme molecule was included and connected by creating a covalent bond between the heme's iron atom and the sulfur atom of conserved cysteine (Cys 432 ) with the propionates of the heme interacting with the side chains of Trp 127 , Arg 131 , and Arg 369 . The CYP397A1V2 model displayed Phenylalanine-cluster (Phe 120 , Phe 121 , Phe 138 and Phe 297 ) above the active site region, with the aromatic side chains stacking against each other to form a prominent hydrophobic core. The Ramachandran plot for CYP397A1V2 showed that approximately 94% of all residues were within the generously allowed region and 1.7% of residues were in the disallowed region (Additional file 13). Thus, the CYP397A1V2 model confirms to the conserved structural folding of P450s, with some unique features including: i) the hydrophobic region was located around the loop (following helix A" and G'-G" helices), ii) the region following helix F appeared shorter compared to other P450s and includes a ordered stretch of polypeptide chain that does not confirm to any secondary structural motif ( Figure 6A), and iii) the catalytic pocket of the CYP397A1V2 displayed a smaller volume compared to other insect P450 proteins ( Figure 6B).

Molecular docking
We used the program Autodock to examine the binding mode of CYP397A1V2 with various insecticides (DDT, imidacloprid, deltamethrin, permethrin and diazinon), any of which C. lectularius may encounter worldwide. The binding constants and free energy change (Additional file 14) revealed DDT, permethrin, deltamethrin and imidacloprid as potential substrates due to their tight-fit into the active site ('V'-shaped hydrophobic cavity) of CYP397A1V2 ( Figure 6C-F). The moiety of DDT and permethrin molecules was located within the binding pocket, and they were adjacent to hydrophobic residues of the Phenylalanine cluster, Val 294 , Ala 298 , Lys 368 , Val 366 and Leu 364 of the derived model CYP397A1V2 ( Figure 7A and 7D). Thus, the interaction of DDT and permethrin with CYP397A1V2 appears to be dominated by mainly hydrophobic interactions. The insecticide substrates deltamethrin and imidacloprid were also fully enclosed by receptor residues in the cavity of CYP397A1V2 and were positioned well within the network of hydrophobic as well as hydrophilic residues. In the case of deltamethrin and imidacloprid, the common receptor residues Tyr 209 , Ala 298 and Thr 302 contributed to the formation of hydrogen bonds between the CYP397A1V2-deltamethrin and CYP397A1V2-imidacloprid complex ( Figure 7B and 7C).
DDT and pyrethroids (permethrin and deltamethrin) are preferential sodium channel modulators in insects [59]. The frequent exposures to DDT, an early generation insecticide during the last century may have altered the development of resistance in C. lectularius to the next generation insecticides through a phenomenon called cross resistance or multi-level resistance [60]. The latter scenario might be the likely phenomenon associated with modern C. lectularius. Our docking effort on CYP397A1V2 supports this hypothesis. Similar observations were found in Papilio glaucus, wherein PcCYP6B4 had the ability to metabolize a broad range of substrates [61]. The unlikely binding of diazinon with CYP397A1V2 suggests potential involvement of other P450s or metabolic proteins that interact with polar compounds.

Validation for kdr mutation
We also examined the PR and PS strains for mutations in the voltage-gated sodium ion channels, as this is a well-characterized kdr mechanism found in C. lectularius [17]. All three PR strains showed both mutations (V419L and L925I) (Additional file 15). Although the PS1 and PS2 strains did not carry any of these mutations, the PS3 strain showed a single mutation at V419L (Additional file 15). These results further suggest that resistance in the PR strains could be attributed to multiple mechanisms of resistance including mutations as well as penetration and metabolic resistance.

Conclusions
In summary, we have developed significant molecular resources for C. lectularius and have identified several candidate genes potentially involved in different phases of insecticide metabolism as well as those participating in other modes of insecticide resistance in this species. Our GLMseq approach further revealed significant differentially expressed ESTs across biologically variant samples. Specifically, the high occurrence of up-regulated cuticular proteins and the expression patterns of P450s in cuticular tissue might represent unique sites for penetration resistance as well as metabolic resistance in C. lectularius. Molecular modeling and docking studies revealed the potential of P450s to metabolize multiple insecticides in C. lectularius. Future functional studies (RNA interference) on cuticular and P450 proteins could lay the foundation for identifying hot-spots for insecticide resistance in C. lectularius, which could provide the basis for developing effective management strategies.

Insect material for RNA-Seq
Six populations (3 pesticide-susceptible strains-Harlan 1, Harlan 2, and FV strain hereinafter referred to as PS1, PS2 and PS3; 3 pesticide-resistant strains collected during 2010 from three residences in different zip codes in Columbus, OH, hereinafter referred to as PR1, PR2 and PR3) of C. lectularius were used in the current study. Bed bug populations were maintained under ambient conditions in the laboratory (~22 ± 1.22°C and 35 ± 6% RH). For each bed bug population, all stages and representatives of multiple generations were housed together on filter paper strips contained in a glass jar (13 cm high by 7 cm diameter) (narrow-mouth Mason pint jar, Ball Corp., Broomfield, CO) with an organza fabric and filter paper covering held in place using a screw-on ring. Multiple jars were needed for large populations. Each bed bug population was fed in situ approximately every 2 weeks on heparinized chicken blood using the Hemotek 5W1 membrane feeding system for blood-sucking insects (Discovery Workshops, Accrington, United Kingdom) with Parafilm ® as the membrane.
The resistance status of bed bug populations was assessed using a discriminating dose adapted from Romero et al. [13]. The Harlan strain was included as a standard susceptible strain given its long-term laboratory rearing status (since 1973). For this test, 10 adult bed bugs from each test population were placed on filter paper discs treated at a rate of 0.13 mg/cm 2 with technical grade deltamethrin (99% purity, Chem Service, Westchester, PA) dissolved in acetone. Acetone-treated filter paper was used as a control. Each population was replicated in triplicate. The condition of the bugs assessed after 24 h exposure to the treated filter paper. At 24 h, 100% of susceptible bed bugs (Harlan strains [PS1, PS2]) were either moribund or dead after exposure to technical deltamethrin, whereas resistant bed bugs (PR1, PR2, PR3) did not show any signs of intoxication. Preliminary resistance testing with PS3 indicated that despite having a known kdr mutation (V419L) (Additional file 15), it was succeptible to dried residues of both lambda-cyhalothrin and deltamethrin at or near the label rate.

RNA isolation and Illumina Paired End (PE) library preparation
Total RNA was extracted separately from each individual (8 adults per population of PS1, PS2, PR1 and PR2) using TRIzol (Invitrogen) and pooled before conducting a quality check (Thermo Scientific Nanodrop 2000 and Agilent Bioanalyzer, Ohio State University, Columbus, OH, USA). Illumina high-throughput sequencing was done using the GAII platform at the Molecular and Cellular Imaging Center, Ohio Agriculture Research and Development Center, Wooster, OH, USA. Samples were prepared for Illumina GA sequencing using the PE library preparation kit (Catalogue Number PE-102-1001, Illumina, San Diego, CA, USA) and TruSeq RNA Sample Prep Kits (Catalogue Number FC-122-1001, Illumina, San Diego, CA, USA) per manufacturer instructions. Briefly, Sera-mag Magnetic Oligo (dT) beads were used for the poly (A) RNA enrichment and divalent cations were used to fragment the purified mRNA (100-400 bp) by heating the reaction mixture at 94°C for 5 minutes. The fragmented RNA was used for double strand synthesis followed by cDNA synthesis. The DNA was end-repaired and phosphorylated using T4 DNA polymerase, Klenow DNA polymerase and T4 PNK. These fragments were 3' adenylated using Klenow Exo-(3' to 5' exo minus) and Illumina PE adapters were ligated using DNA Ligase. To select appropriate size (for PE library preparation kit) and to eliminate unligated adapters, the adapter ligated products were puri-

Raw data analysis and de novo assembly
To eliminate low quality nucleotides, raw Illumina reads were trimmed using a custom Perl script with windowed adaptive trimming (Phred quality threshold of 20 and minimum read length of 20 nt). To do the de novo assembly, the trimmed reads from PR sample 1 and PS sample 1 were separately fed into Rnnotator, automated de novo RNAseq assembly pipeline [62], to remove duplicates and erroneous sequences. Multiple rounds on velvet assemblies then were performed using different hash values to account for the different sequencing depth for different transcripts [62]. Resulting contigs were merged using Minumus2 from the AMOS package [62,63]. These were combined with previously published C. lectularius 454 sequences [19] for better coverage and quality using Minumus2 from the AMOS package to make the final reference contigs.

Validation and annotation of assembled contigs
To validate the unigenes, we used three criteria: BLASTx searches (E value < 10 -3 ) between unigenes and the NCBI non-redundant protein (nr) database http://www.ncbi. nlm.nih.gov; BLASTx searches (E value < 10 -3 ) between unigenes and several insect species with genome sequence information (A. pisum, A. gambiae, D. melanogaster, P. humanus); and prediction of putative open reading frames (ORFs) using http://www.scbi.uma.es/cgibin/full-lengther/full-lengther_login.cgi. The ortholog hit ratio was calculated using custom python script [64]. To annotate the unigenes, BLASTx searches (E value < 10 -3 ) were performed between nr and unigenes. BLAST results were imported to Blast2GO program for further annotation of the unigenes [65]. After the mapping step, those gene ontology (GO) terms with E value < 1E-3, annotation cut-off > 45, and GO weight > 5 were used for annotation. To find the pathways in which putative peptides of the unigenes are involved, analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed using Blast2GO [66]. To find enriched GO terms, enrichment analysis (Fisher's Exact Test) tool in Blast2GO software was used with term filter value below 0.05, term filter mode "FDR" and two-tailed test. To categorize the GO terms into different GO categories, a web-based tool, CateGOrizer with "Aqua" tool was used and these categories were compared to a precomputed GO terms from Drosophila melanogaster [67,68].

Differential gene expression analyses
Differential gene expression analysis among PR and PS samples was performed using two-way GLMseq, a custom R-script for fitting a two-way Generalized Linear Model to the four un-normalized RNA-Seq counts at each contig. The GLM function in R was used to estimate and test both the replicate effect and the resistance effect for each quartet of Poisson counts, with offsets determined by the total mapped reads for each sample. This model-free method makes no assumptions about the four unknown Poisson parameters and requires no pseudo-reference for the calculation of p-values.
The model treated each count as an independently measured Poisson variate with offset determined by library size. The examples below illustrate the effectiveness of GLMseq for setting aside secondary batch effects: In example A (Contig_18664), the raw counts for Batch-1 (4, 4) and Batch-2 (75359, 22) clearly suggest a large interaction effect (Additional file 16). In this case the GLMseq estimate of the main Resistant/Susceptible fold-change was 1.14, with a two-sided p-value of 0.85, indicating the absence of a significant resistance effect. In example B (Contig_27586), the raw counts for Batch-1 (18,4) and Batch-2 (15041, 527) clearly suggest a positive resistance effect (Additional file 16). In this case the GLMseq estimate of the main Resistant/Susceptible fold-change was 5.15, with a two-sided p-value of 0.003, indicating the presence of a significant positive resistance effect. In both examples, naïve analysis of the pooled Resistant and Susceptible counts would have mistakenly assigned much lower p-values to the estimated fold-changes.
For finding uniquely expressed genes in PR and PS, the read counts were normalized using EDASeq method, which performs normalization for sequence length and GC content for intra and inter samples [69]. We considered ESTs with less than 10 short reads aligned in both biological replicates as not reliably expressed; for the expressed genes, we used more than 100 short reads aligned in both biological replicates as an arbitrary cut off.

cDNA preparation for quantitative PCR analysis
Individuals from the six strains (3 PR and 3 PS) of C. lectularius were categorized into various development stages (early instar nymphs, late instar nymphs and adults) as previously described [19,70]. The tissues (cuticle, midgut and Malphigian tubules) were dissected as per Mamidala et al. [71]. Total RNA was extracted using TriZOL and evaluated with Nanodrop. Further, RNA was treated with TURBO DNase™ (AMBION, Inc., Austin, TX) to remove any genomic contamination, and stored (at -80°C) until further use. First strand cDNA was synthesized in 20 μL reactions using~0.5 μg of RNA, oligo dT primer and Super Script III First-Strand Synthesis Super Mix (Invitrogen) and the resultant cDNA was diluted to 20 ng/μl for qRT-PCR studies.

Gene characterization of CYP 397A1V2
To develop full length CYP397A1V2, we performed 3' RACE (Rapid Amplification of cDNA Ends) using an oligo dT primer (Invitrogen) and gene specific primers (Additional file 17). The amplicon was sent for sequencing at the Molecular and Cellular Imaging Center, Wooster, OH, USA. Two additional gene specific primers were designed to the ends of the full length gene to confirm the sequence (Additional file 17). The identification and annotation of the CYP3 clan member was performed using BLASTx search against nr database at NCBI http://www.ncbi.nlm.nih.gov/. For the phylogenetic analysis, an unrooted neighbor-joining tree was constructed with 500 bootstrap replicates and excluding positions with gaps using MEGA version 5 [72].

Primer designing and qRT-PCR
Twelve potential genes involved in insecticide resistance that showed differential expression in RNA-Seq were profiled including 4 cytochrome P450s (CYP9, CYP6A2, CYP6A13 and CYP301A2), 5 cuticular proteins (larval cuticle protein, pupal cuticle protein, cuticle protein analogous to peritrophin, chitin deacetylase and chitin synthase), an antioxidant gene (superoxide dismutase), an acetylcholinesterase gene, and an ABC transporter. Primers were designed using Beacon Designer 7 software (Additional file 17). Prior to qRT-PCR, standard PCR was performed for all primer pairs and their products were run on agarose gel electrophoresis to ensure single bands. qRT-PCR reactions were performed in 96well plates using a BioRad thermocycler (CFX-96) as per Bai et al. [19]. The analysis included three biological replicates and three technical replicates (within each biological replicate) for RNA-Seq validation and for developmental stages whereas two biological replicates (PS1, PS2, PR1 and PR2) and two technical replicates (within each biological replicate) were included for tissues. Relative gene expression among PR and PS adult populations was analyzed as per Bai et al. [19]. For tissues and developmental stages, the fold change in gene expression between PR and PS strains of C. lectularius were derived by the 2 -ΔΔCT method [73] using ribosomal protein (RPL18) of C. lectularius as the internal control gene [71].

Statistical analysis
The Wilcoxon rank sum test revealed statistically significant differences (P < 0.0001) between the adult samples for PR and PS strains [74]. The transformed C T values (2 -C T ) were used for statistical analysis using SAS (SAS/ STAT User's Guide, Version 9.1, SAS Institute Inc.) with a significance level (α) of 0.05. A PROC MIXED analysis of variance (ANOVA) was performed for tissues and developmental stages of PR and PS strains.

Construction of the CYP397A1V2 model
The MODELLER 9v8 program [75] was used to construct the CYPP397A1V2 structure. MODELLER is a general program that implements comparative protein structure modeling by satisfying spatial restraints in terms of probability density functions. To this end, the crystal structure of Human Cytochrome P450 3A4 (PDB code 1W0E) was used as a 3D template (Additional file 12). The MODELLER program was applied to generate 40 satisfactory models, including the ironoxo group of CYP397A1V2. The model with the lowest energy and the lowest restraint violation was selected. The initial model was improved by energy minimization. First, energy minimization of 1,000 steps of steepest descent (SD) followed by 2,000 steps of conjugate gradient (CG) was carried out in order to release the conflicting contacts among residues. After performing 2,000 steps of conjugate gradient (CG) minimization, MD simulation was carried out to examine the quality of the model structure, by checking the stability via performing 2000 ps simulations at a constant temperature 300 K. An explicit solvent model SPC216 water was used [76]. All energy minimization and MD simulations were accomplished by GROMACS4.0 software package [77,78] using the GROMOS96 43a1 force field [79,80].

Molecular docking, preparation of the protein and the ligand
The molecular docking program Autodock, which uses a powerful Lamarkian genetic algorithm (LGA), was used to dock the ligands to the protein active site [81]. Insecticide substrate models (DDT, diazinon, deltamethrin, imidacloprid and permethrin) were built, and their geometry was optimized through discover3 in the InsightII/ Builder program http://www.accelrys.com. To recognize the binding sites in CYPP397A1V2, blind docking was carried out with the grid size set to 126, 126 and 126 along the X, Y and Z axes with 0.397 Å grid spacing. The AutoDocking parameters used were GA population size: 150 and maximum number of energy evolutions: 250,000. During docking, a maximum number of 10 conformers was considered, and the root-mean-square (rms) cluster tolerance set to 1.5 Å. One of the lowest energy conformations was considered for further analysis.

DNA extraction and gene sequencing for kdr mutations
For each of the bed bug populations (PS1, PS2, PS3, PR1, PR2, PR3), five whole adults were homogenized and genetic material was extracted using the EZNA DNA extraction kit (Omega Bio-Tek). Gene segments containing known mutation sites in the voltage-gated sodium ion channel were amplified using PCR as described in Zhu et al. [17] (Additional file 17). PCR products were subsequently purified and sequenced by Functional Biosciences (Madison, WI, USA) http://www. functionalbio.com/web/. with a two-sided p-value of 0.85, indicating the absence of a significant resistance effect, whereas the Contig_27586 clearly suggests a positive resistance effect.
Additional file 17: Primers used in the current study. List of oligonucleotide primers used in qRT-PCR validation.