CITRUS RED MITE REARING AND RESISTANT MITE SELECTION

A colony of citrus red mite, P. citri, was collected from citrange (Citrus sinensis × Poncirus trifoliata ) next to citrus orchard fences (Citrus Research Institute, Chinese Academy of Agricultural Sciences), which had not been exposed to acaricides for more than 10 years. These mites were reared on citrange for 3 years under acaricide-free conditions (25 ± 1°, 80 ± 5% RH and 14:10 h L: D and were considered to be a susceptible strain (SS). The SS strain was treated with hexyzhiazox (about a 70% population mortality rate was observed at every spray application) and was continuously screened for 20 generations until it was regarded as a fully resistant (RR). The selection process and the resulting changes in resistance were shown in a previous study (Liu et al. 2011a). A mixture of eggs, larvae and adults was collected in PerkinElmer (PE) tubes with one sample for SS and RR respectively. The 2 samples were stored at -80° until use.

RNA ISOLATION AND INTEGRITY EXAMINATION

Total RNA was isolated from the following 2 samples: RR and SS of P. citri to hexythiazox. From each sample, approximately 8 mg of mites were homogenized with liquid nitrogen in a mortar. RNA was extracted with TRIzol® reagent (Invitrogen, USA) according to the manufacturer's instructions and was treated with RNase-free DNase I (Takara Biotechnology, China). RNA integrity was confirmed with a minimum RNA integrated number value of 8 by the 2100 bioanalyzer (Agilent).

DIGITAL GENE EXPRESSION (DGE) LIBRARY CONSTRUCTION AND SEQUENCING

Approximately 10 μg of total RNA from each specimen (a mixture of RNA from eggs, third-instar larvae, pupae, and adults at equal ratios) was used to construct the DGE libraries. Poly(A) mRNA was isolated with oligo-dT beads and then treated with the fragmentation buffer. The cleaved RNA fragments were then transcribed into first-strand cDNA using reverse transcriptase and random hexamer-primers. This was followed by second strand cDNA synthesis using DNA polymerase I and RNase-H. The double stranded cDNA was further subjected to end-repair using T4 DNA polymerase, Klenow fragment, and T4 polynucleotide kinase followed by a single <A> base addition using Klenow 3′ to 5′ exo-polymerase. It was then ligated with an adapter or index adapter using T4 quick DNA ligase. Adaptor ligated fragments were selected according to size and the desired range of cDNA fragments were excised from the gel. We performed a PCR reaction to selectively enrich and amplify the amount of the fragments. Finally, after validation on an Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System, the DGE libraries was sequenced on a flow cell using an Illumina HiSeq2000.Bioinformatics analysis process of sequence data were presented in Supplemental File 1. In particular, the raw reads were cleaned up by removing adapter sequences, reads in which unknown bases are more than 10% and low quality reads of quality value ≤ 5 as described by Liu (Liu et al. 2011a).

MANUAL CURATION OF CYPS IN P. CITRI

We blasted P. citri transcriptomes with sequences from the NCBI nucleotide database and obtained genes related to general pesticide resistance in other insect species. In total, there were 121 related to cytochrome P450 monooxygenase (CYPs). We chose only 46 of these for functional annotation according to sequence matching (E-value< le-5) owing to low sequence similarity of other P450s. Nucleotide sequences were dynamically translated using the EXPASY Proteomics Server ( http://www.expasy.ch/tools/dna.html, Swiss Institute of Bioinformatics). All the identified sequences were searched with BLASTx against all the assembled contigs in the iceblast server using an E-value cut-off of 1e-5 and the results with more than 99% similarity with the query sequence were eliminated as allelic variants (note that from those sequences, only the longest contigs with the best coverage were manually curated). We also searched the P. citri P450s on the NCBI network and obtained 2 amino acid alignments of P. citri CYPs. Next, we manually curated each of the CYPs with the help of comparison with the CYP genome of the two-spotted spider mite ( http://drnelson.uthsc.edu/Two.spotted.spider.mite.htm). Then, we assigned gene names based on the homology of P. citri CYPs to T. urticae using defined nomenclature and naming rules (Nelson 2006).

PHYLOGENETIC ANALYSIS OF CYPS IN P. CITRI

MEGA 4.0 software (Tamura et al. 2007) was used to analyze the phylogenetic relationships between P450 of P. citri and several CYPs in T. urticae to predict their classification. The neighbor-joining method was used to create phylogenetic trees (Saitou & Nei 1987). Positions containing alignment gaps and missing data were eliminated with pairwise deletion. Bootstrap analysis of 1,000 replication trees was performed to evaluate the branch strength of each tree (Efron et al. 1996). The NJ tree included bootstrap confidence levels at nodes, which reflected the confidence of trees sampled during the reconstruction that included each particular branch. Only clades with a bootstrap value higher than 50 were selected for the bootstrap consensus tree (Soltis & Soltis 2003).

SEQUENCE ALIGNMENT AND VERIFICATION OF MUTATIONS IN CYPS IN P. CITRI

Multiple nucleotide sequence alignment of each CYP between the resistant and susceptible haplotypes of P. citri was performed using ClustalX (1.83) software (Thompson et al. 1994) and followed methods of Liu et al. 2011a. This alignment process was also based on the Feng & Doolittle (1987) algorithm, but with an improved choice of alignment parameters using dynamic assignment of penalties. Two pairs of specific Primers (Table 2) for surveying the nucleotide differences between RR and SS were designed by using Primer Premier 5.0 software based on the transcriptome of P. citri. Analyses using PCR were carried out in a total volume of 25 μL and performed in a thermocycler(GeneAmp PCR system 9600; Perkin-Elmer Corporation, Norwalk, Connecticut). The PCR reaction contained 12.5 μL of 1×Taq MasterMix from (Tiangen Biotech), 10 pmol of each primer, 30 ng of cDNA. PCR was performed using a Mastercycler® Gradient using the following protocol: initial denaturation at 94 °C for 3 min, followed by 30 cycles of 94 °C for 30 s, 56 °C for 30 s, 72 °C for 30 s, with a final extension step of 72 °C for 10 min. Samples (5 μL) of reaction mixtures were examined by electrophoresis through 2% agarose gels in TAE buffer. Bands were revealed by visualization with UV light after ethidium bromide staining. PCR products were purified and sequenced by Shanghai Sangong Biological Engineering Technology & Service CO., LID. All these primers mentioned above were also used as sequencing primers. Once a mutation was confirmed, both nucleotide sequences from the 2 strains were dynamically translated using the EXPASY Proteomics Server ( http://www.expasy.ch/tools/dna.html, Swiss Institute of Bioinformatics). Amino acid sequence alignment was performed by ClustalX (1.83) software to confirm whether the mutation was a sense mutation.

DIFFERENTIAL EXPRESSION OF CYPS IN P. CITRI

The uniquely mapped reads for a specific transcript were counted by mapping reads to assembled sequences using SOAP (Li et al. 2009). Then the RPKM value for each transcript was measured in reads per kilobase of transcript sequence per million mapped reads (Mortazavi et al. 2008). The transcript fold change was calculated by the formula of log₂ (RR_ RPKM / SS_RPKM). If the value of either RR_ RPKM or SS_RPKM was zero, we used 0.01 instead of 0 to calculate the fold change. We modified Audic's (Audic & Claverie 1997) method to analyze differential expression. The probability of a specific gene being expressed equally between the 2 samples was defined by the following formula:

In this equation, N1 and N2 indicate the total number of clean reads in the SS and RR haplotypes and x and y represent the mapped clean read counts of one transcript in the 2 samples. The false discovery rate (FDR) method was used to determine the threshold of the P value in multiple tests. In this study, we used ‘FDR ≤ 0.001 and the absolute value of the log₂ ratio ≥ 1’ as the threshold to judge the significance of differential gene expression.

QUANTITATIVE RT-PCR ANALYSIS

Quantitative Reverse transcription polymerase chain reaction (RTPCR) detection was used to verify 5 differentially expressed genes (2 up-regulated and 3down-regulated genes), e.g., CYP389A6, CYP385A2, CYP307A1, CYP307A2, CYP389A1. Multiple specific Primers (Table 2) were designed by using Primer Premier 5.0 software (PREMIER Biosoft International, California). The expression study was performed on an Mx3000P™ Multiple Quantitative PCR System (Stratagene, La Jolla, California, USA) with ELF1A, elongation factor-1 alpha, used as reference gene (Niu et al. 2012). The reaction system and reaction parameters were same as Niu et al. (2012) described earlier. After collecting the Ct values, the relative quantitative expression of these 5 genes were calculated by the 2^-δδCT method (Livak & Schmittgen 2001).

DIGITAL GENE EXPRESSION (DGE) LIBRARY SEQUENCING

Based on the transcriptome sequence data, 2 DGE libraries were constructed to identify the expression profiles of the various strains. After removing the low-quality reads, each library generated more than 7 million clean reads. Among these clean reads, approximately 2.8 million and 2.9 million (36.92% and 40.36%, respectively) were mapped to unigenes in each library (Table 3). The percentage of clean reads ranged from 98.68% to 99.70%, reflecting a high quality of sequencing.

MANUAL ANNOTATION IDENTIFIED 48 CYPS IN P. CITRI

Manual annotation and curation of the CYPs in the P. citri transcriptome sequence assembly ( http://www.ebi.ac.uk/ena/data/view/ERP000885) produced 121 sequences related to cytochrome P450 monooxygenase (CYPs). Of these, 46 were manually curated, as the remainder was found to have either low sequence similarities to other P450s according to sequence matching (E-value< le-5), or they contained too many sequencing errors. In addition, 2 full length P450s were detected by Jiang et al. (Jiang et al. 2010), which were the first detected P450 genes of P. citri. The latter were added to our data, and these 48 P450 sequences were named by D. R. Nelson in accordance with the P450 nomenclature committee conventions ( http://drnelson.uthsc.edu/cytochromeP450.html). Based on the closest BLASTX matches in the NCBI nr database first, and when possible, by phylogenetic analyses with CYPs in T. urticae ( http://drnelson.uthsc.edu/Two.spotted.spider.mite.htm), the P450s of P. citri were assigned to all 4 major insect CYP clans (CYP2, CYP3, CYP4 and mitochondrial). The 4 clans can be further subdivided into 15 families, 24 subfamilies, and 48 putatively functional isoforms (Fig. 1). Comparison of the number of functional CYP genes in different transcriptomes and genomes (Table 1) shows that the CYP number is lower than in the typical invertebrate or insect with the exception of the honeybee (Claudianos et al. 2006) and the New World screwworm (Carvalho et al. 2010).

Fig. 1.

Number, family and clan distribution of cytochrome P450 genes in Panonychus citri. The number shown along each column represents the P450 family and the number in parenthesis is the number of individual genes in the corresponding family. The P450 gene sequence information generated is from the VectorBase of the P. citri transcriptome sequence.

Table 1.

Comparison of the number of functional CYP genes in different genomes or transcriptomes of various arthropods.

PHYLOGENETIC ANALYSIS OF CYPS IN P. CITRI.

Repeated and exhaustive searching for P. citri transcriptome assemblies and nucleotides in the NCBI uncovered a total of 48 CYP genes. Based on inferred amino acid sequences, these genes fell into 14 CYP gene families. Molecular phylogenetic analysis (Fig. 2) by the neighborjoining method (Saitou & Nei 1987) shows the relationships among CYP genes and gene families in P. citri and several T. urticae P450s. Positions containing alignment gaps and missing data were eliminated with pairwise deletion. Bootstrap analysis of 1,000 replication trees was performed to evaluate the branch strength of each tree (Efron et al. 1996). The NJ tree included a bootstrap confidence level at nodes that reflects the confidence of trees sampled during the reconstruction that included each particular branch. In general, the tree demonstrates that the 4 major clans found in insects (Feyereisen 2006), which include CYP2, CYP3, CYP4 and mitochondrial clans, encompass all of the CYPs in P. citri. All nucleotide and amino acid sequences of cytochrome P450 monooxygenase genes in P. citri are available in the database ( http://www.ebi.ac.uk/ena/data/search?query=Panonychus citri).

Fig. 2.

Neighbor-joining phylogenetic analysis of cytochrome P450 from Panonychus citri and Tetranychus urticae. 4clans were observed. There are species (P. citri and T. urticae) in the phylogenetic tree. Only 10 sequences belong to T. urticae; A (Pc) before the CYP name denotes P. citri, a (Tu) before the CYP name denotes T. urticae. Numbers at nodes are bootstrap values.

Table 2.

Specific RNA primers used in sequencing and quantitative RT-PCR of Panonychus citri.

Fig. 3.

Quantitative Real-time PCR analysis of CYPs in Panonychus citri between the hexythiazox-resistant (RR) and susceptible (SS) strains. The numbers of genes down-regulated and up-regulated in the RR relative to the SS are indicated above or below the X axis. The light or dark gray was susceptible strain and resistant strain, respectively.

The CYP2 clan of P. citri contains 12 members, which can be separated into 2 distinct families (CYP307, CYP392). The spook gene CYP307A1, which is expressed in the prothoracic gland (Niwa et al. 2005), is a conserved Halloween gene involved in the early stages of ecdysone synthesis (Rewitz et al. 2006; Rewitz & Gilbert 2008), and can convert 7-dehydrocholesterol to δ4-diketol with the gene product of spookier (CYP307A2) (Gilbert 2008). These 2 genes are close paralogs that are believed to mediate the same enzymatic reaction at different stages of development (Rewitz & Gilbert 2008). The other CYP2 clan family in P. citri (CYP392) is divided into 3 subfamilies.

Fig. 4.

Nucleotide sequence comparison of the CYP307A1 in Panonychus citri between the hexythiazox-resistant (RR) and susceptible (SS) strains. Mazarine shading indicates identities and different color shading represents mutations. “-” represents no sequence to compare. Three SNP sites were detected in all. The first nucleotide mutation (A to C) is located at 841, the second mutation is 1395-T to C, and the final mutation is 1491-T to C. . This figure is shown in color in a supplementary document online as Suppl. Fig. 4 in Florida Entomologist 98(1) (March 2015) at  http://purl.fcla.edu/fcla/entomologist/browse.

The CYP3 clan of P. citri contains 12 CYPs. Genes in the CYP3 clan are the most numerous among insect P450 genes, and are often found in large clusters associated with oxidative detoxification of xenobiotics (Kretschmer & Baldwin 2005; Li et al. 2009; Waxman 1999) and endobiotics (Strode et al. 2008; Zhang et al. 2004). The first insect P450 gene (CYP6A1) was isolated from an insecticide-resistant strain of the housefly, M. domestica L. (Feyereisen et al. 1989). The first P450 gene (CYP6A2) was isolated from D. melanogaster, and is also associated with resistance to insecticides (Waters et al. 1992). In P. citri, the CYP3 clan contains 12 genes and 5 subfamilies. Eight genes belong to the CYP385 family and just 1 gene (CYP382A1) is part of the CYP382 family.

Table 3.

Alignment statistics of the DGE-SEQ analysis of Panonychus citri RNA extracted from the hexythiazox-resistant (RR) and susceptible (SS) strains.

The CYP4 clan of P. citri contains 19 genes and can be divided into 6 families and 10 subfamilies. Interestingly, the number of CYP4 clan genes varies in different species. For example, there are 32 in D. melanogaster (Tijet et al. 2001), 45 in Anopheles gambiae Giles (Holt et al. 2002), 44 in Tribolium castaneum (Herbst) (Richards et al. 2008), 32 in B. mori (L.) (Bin et al. 2005), 59 in Aedes aegypti L. (Strode et al. 2008), and only 4 in Apis mellifera L. (Claudianos et al. 2006). The great diversity of genes in the CYP4 clan was also reflected in a great diversity of functions. CYP4s in other insects have been implicated in functions as diverse as 20-hydroxyecdysone biosynthesis (Maibeche-Coisne et al. 2000), pheromone metabolism (Maibeche-Coisne et al. 2000) and pyrethroid insecticide resistance (Pridgeon et al. 2003).

The mitochondrial clan of P. citri probably originated in the CYP2 clan, which includes P450s with essential physiological functions, and contains 5 members in 3 families and 3 subfamilies. Two of the members, CYP302A1 and CYP314A1, are highly conserved halloween genes, involved in ecdysone synthesis (Rewitz et al. 2006, 2008). The disembodied (dib) gene (CYP302A1) codes a cytochrome P450 enzyme that adds a hydroxyl group to the carbon-22 position of 2, 22, dE- Ketotriol to make 2- Deoxyecdysone (Gilbert 2004). The dib mutants are defective in producing cuticle and have severe defects in morphological processes such as head involution, dorsal closure and gut development (Chávez et al. 2000). The shade gene (CYP314A1) codes for a P450 enzyme which reportedly adds hydroxyl group to the 2C- position of ecdysone to make 20-hydroxyecdysone, as the final step in the biosynthetic pathway (Petryk et al. 2003). CYP381 was a single family in P. citri.

Fig. 5.

Alignment of the predicted amino acid sequences of CYP307A1 in Panonychus citribetween the hexythiazox-resistant (RR) and susceptible (SS) strains. Mazarine shading indicates identities and different color shading represents mutations. “-” represents no sequence to compare. Only one amino acid mutation (278-lysine to glutamine) was detected. . This figure is shown in color in a supplementary document online as Suppl. Fig. 5 in Florida Entomologist 98(1) (March 2015) at  http://purl.fcla.edu/fcla/entomologist/browse.

SEQUENCE ALIGNMENT AND VERIFICATION OF MUTATIONS IN CYPS OF P. CITRI

ClustalX (1.83) software (Thompson et al. 1994) was used to analyze alignment differences in each CYP between the RR and SS of P. citri. After sequence verification, 3 SNP sites were detected inCYP307A1 (1688 bp, Fig. 4). The first nucleotide change located at 841(A to C), the second change located at 1395 (T to C), and the final change located at 1491 (T to C). With alignment of the 2 predicted amino acid sequences of the CYP307A1 (Fig. 5), a mutation was detected for one amino acid (278-lysine to glutamine). Therefore, the mutation was a sense mutation and the other 2 were nonsense mutations. A single SNP site, 40 (A to T), was detected in CYP381A2 (619 bp, Fig. 6). With alignment of the 2 predicted amino acid sequences of the CYP381A2 (Fig. 7), there was a sense amino acid mutation detected (14-threonine to serine). In insects, steroidogenic CYPs are products of the Halloween genes phantom (CYP306A1), disembodied (CYP302A1), shadow (CYP315A1) and shade (CYP314A1) and are responsible for the last 4hydroxylations in the pathway leading to 20E (Niwa et al. 2004; Niwa et al. 2005; Warren et al. 2002),which are biochemically similar to one that yields 20E in crustaceans (Lachaise et al. 1993). In D. melanogaster, mutations in these genes disrupt 20E production and cause the arrest of embryonic development and death. Spook (CYP307A1) is another member of this CYP group and, when mutated, results in low 20E mutants (Niwa et al. 2005). In B. mori, the nucleotide transition of C to T at position 1049 was found in phantom (CYP306A1), resulting in a nonsense mutation at amino acid position 286. The morphology of the phantom mutant is normal until stage 14 when there is a failure of dorsal closure and head involution and embryos become compacted. The phantom mutant exhibited severe reductions in the epidermal expression of these 20- hydroxyecdysone-inducible genes (Niwa et al. 2004).

DIFFERENTIAL EXPRESSION OF CYPS IN RR AND SS OF P. CITRI

A rigorous algorithm (see Methods) was developed to identify genes differentially expressed between resistant (RR) and susceptible (SS) mite populations treated with hexyzhiazox by referring to “The significance of digital gene expression profiles” (Audic & Claverie 1997). To explore the gene expression levels further, the reads per kb per million reads (RPKM) method (Mortazavi et al. 2008) was adapted to eliminate the influence of variation in gene length and the total reads number. Results showed 22 P450s were up-regulated and 24 were down-regulated (Table 4). Furthermore, the significance of gene expression differences were determined using the false discovery rate (FDR≤ 0.001) and the absolute value of the log₂ ratio (≥ 1). Nine genes were differentially expressed between the 2 samples at significant levels, including 5 down-regulated and 4 up-regulated genes (Table 4). Seven members of P450s had a greater value than 2 based on a log₂ ratio formula, and these members were distributed as 2, 2 and 3 in the CYP 2 clan, CYP 3 clan and CYP 4 clan, respectively. It was straightforward to identify the highest up-regulated expression gene (log₂ ratio = 12.14) as CYP389A6, which belonged to the CYP4 clan. The gene of highest expression might be analogous with those CYP4s in other insects that have been implicated in pesticide metabolism. For example, CYP4G8 is over expressed in pyrethroid-resistant strains of H. armigera (Pittendrigh et al. 1997). CYP4C27 in A. gambiae is over expressed in a DDT-resistant strain (David et al. 2005) and CYP4G19 expression in Blattella germanica is correlated with pyrethroid resistance (Pridgeon et al. 2003). Several CYP4 genes are over expressed in pesticide resistant C. pipiens and Diabrotica virgifera (Scharf et al. 2001; Shen et al. 2003). The expression level of CYP385A2 (log₂ ratio = 9.46), belonging to the CYP3 clan, was a little lower than that of CYP389A6. The CYP3 clan is related to metabolism of pyrethroids, chlorinated hydrocarbons (e.g., DDT) and neonicotinoids (e.g., imidacloprid) (Daborn et al. 2002; Jiang et al. 2010; Karunker et al. 2008; Komagata et al. 2010; Nikou et al. 2003). The most down-regulated gene (log₂ ratio = -11.816) was CYP389A1, which belonged to the CYP4 clan. The second lowest (log₂ ratio = -11.079) was CYP307A2, which belonged to the CYP2 clan. The genes with the absolute value of the log₂ ratio (≥ 2) might be related to the lower activity of hexythiazox-poisoned P. citri individuals (Liu et al. 2011a). However, quantitative RT-PCR results showed that the upregulated value of log₂ ratio (RR/SS) in CYP385A2 and CYP307A1 were only 0.28 and 0.17, the down-regulated value of log₂ ratio (RR/SS) in CYP389A6, CYP307A2 and CYP389A1 were just -0.69, -1.46 and -1.08, respectively (Fig. 3). All sequences at  http://www.ebi.ac.uk/ena/data/search?query=Panonychus citri

Fig. 6.

Nucleotide sequence comparison of CYP381A2 in Panonychus citri between the hexythiazox-resistant (RR) and susceptible (SS) strains. Mazarine shading indicates identities and different color shading represents mutations. “-” represents no sequence to compare. Just one SNP site was detected. The nucleotide transition of A to T was at position 40. . This figure is shown in color in a supplementary document online as Suppl. Fig. 6 in Florida Entomologist 98(1) (March 2015) at  http://purl.fcla.edu/fcla/entomologist/browse.

Fig. 7.

Alignment of the predicted amino acid sequences of the CYP307A1 in Panonychus citri between the hexythiazox-resistant (RR) and susceptible (SS) strains. Mazarine shading indicates identities and different color shading represents mutations. “-” represents no sequence to compare. We detected a sense amino acid mutation (14-threonine to serine). . This figure is shown in color in a supplementary document online as Suppl. Fig. 7 in Florida Entomologist 98(1) (March 2015) at  http://purl.fcla.edu/fcla/entomologist/browse.