Introduction
The citrus red mite, Panonychus citri (McGregor), is the most common pest in citrus-producing regions. Panonychus citri has piercing and sucking mouthpart to suck the juice from citrus leaves, tender shoots, flower buds, and fruits, causing defoliation and abscission. These mites spin webs on the citrus to immobilize themselves, preventing the mites from being washed away by rainwater. Moreover, the use of a web expands the feeding range, with the help of wind. Female adult mites produce haploid males with no mating to maintain the gender ratio of the population. The mite population expands rapidly in the spring and autumn, negatively impacting citrus production.
The use of chemical acaricides is the primary strategy for pest management. The Insecticide Resistance Action Committee (IRAC) divides chemical acaricides into 14 categories (Van Leeuwen et al., Reference Van Leeuwen, Tirry, Yamamoto, Nauen and Dermauw2015) depending on the mode of action. Misuse of any category of acaricides may lead to resistance in mites. Panonychus citri has developed varying degrees of resistance to many acaricides, including spirodiclofen, bifenazate, hexythiazox, pyridaben, abamectin, fenpropathrin, and azocyclotin (Yamamoto et al., Reference Yamamoto, Yoneda, Hatano and Asada1996; Ran et al., Reference Ran, Chen and Wang2009; Hu et al., Reference Hu, Wang, Wang, You and Chen2010; Van Leeuwen et al., Reference Van Leeuwen, Van Nieuwenhuyse, Vanholme, Dermauw, Nauen and Tirry2011; Ouyang et al., Reference Ouyang, Montez, Liu and Grafton-Cardwell2012). According to IRAC's classification, amitraz is an octopamine receptor agonist, which causes persistent excitement of insects until death. Unlike low-toxicity acaricides such as spirocyclic tetronic acids (spirodiclofen) and mite growth inhibitors (etoxazole, clofentezine, hexythiazox, and flufenoxuron), amitraz is a wide-spectrum insecticide, which was once rarely used in orchards.
The resistance problem in P. citri has been increasing in recent years. Few new acaricides have been developed. Thus, amitraz has been increasingly applied in orchards to solve this problem, as an alternative to low-toxicity acaricides. However, P. citri is already developing resistance to amitraz (Ran et al., Reference Ran, Chen and Wang2009). Only a limited number of studies have been devoted to the mechanism of the amitraz resistance in mites, as opposed to a great many studies in other parasites. Through calculating gene mutation frequencies, Corley et al. found that the I61F mutation in the β-adrenergic octopamine receptor gene of Rhipicephalus microplus correlates with amitraz resistance (Corley et al., Reference Corley, Jonsson, Piper, Cutullé, Stear and Seddon2013). Based on differential expression analysis in cell lines, Koh-Tan et al. found that the β-adrenergic octopamine receptor and ATP-binding cassette (ABC) B10 might also be involved in the formation of amitraz resistance in R. microplus (Koh-Tan et al., Reference Koh-Tan, Strachan, Cooper, Bell-Sakyi and Jonsson2016). Through eukaryotic protein expression and metabolism, Kita et al. found that DPMF, the metabolite of amitraz, was more likely to activate α- and β-adrenergic-like octopamine receptors (Kita et al., Reference Kita, Hayashi, Ohtani, Takao, Takasu, Liu, Ohta, Ozoe and Ozoe2017). Multiple mechanisms might be involved in amitraz resistance, including base mutations in the β-adrenergic octopamine and the octopamine/tyramine receptors and altered activities of monoamine oxidases and ABC transporters (Jonsson et al., Reference Jonsson, Klafke, Corley, Tidwell, Berry and Caline2018). The results of differentially expressed genes (DEGs) and modeled KEGG pathway in amitraz-treated Rhipicephalus decoloratus laid the basis for further investigation on resistance mechanisms (Baron et al., Reference Baron, Barrero, Black, Bellgard, van Dalen and Maritz-Olivier2018). The stability of cell line from Ixodes ricinus was affected after treatment with amitraz, and ABC transports of cell line were upregulated or downregulated to varying degrees, indicating that this cell line was fit for in vitro mechanistic studies of amitraz resistance (Mangia et al., Reference Mangia, Vismarra, Genchi, Epis, Bandi, Grandi, Bell-Sakyi, Otranto, Passeri and Kramer2018). Functional analyses through in vitro expression of β-adrenergic-like octopamine receptor in Plutella xylostella, Huang et al. found that amitraz, octopamine, and tyramine all activate the β-adrenergic-like octopamine receptor, leading to elevated cyclic AMP levels. An expression pattern study indicated that the β-adrenergic-like octopamine receptor is involved in nearly every growth stage of P. xylostella (Huang et al., Reference Huang, Ma, Deng, Zhu, Liu, Zhou and Zhou2018).
The amitraz resistance mechanism in P. citri remains unclear. We would like to know which proteins or genes play roles in metabolic detoxification of P. citri against amitraz. Once the background of resistant strains is obtained, genetic analysis of resistance can be conducted through forward genetics to attain preliminary date about resistant trait in P. citri on the macroscopic scale. This date will be the basis for localization and analysis of the target genes in the population. Furthermore, studying synergistic effects is helpful in screening for genes related to metabolic resistance and for narrowing the scope of detoxification-related genes. From DEG analysis, we can screen candidate genes consistent with the synergism studies at the molecular level. The next step is to perform in vitro metabolism analyses and genome re-sequencing, to better understand resistance mechanisms.
Materials and methods
Indoor screening for resistant strains
The acaricides-sensitive strain (SS) of P. citri used in the present study was isolated from a wild citrange collected near the orchard belonging to the Citrus Research Institute (Chongqing, China) in 2005. The population was transferred to a pest-free fructus citri aurantii and bred in the laboratory without pesticide exposure for 13 years. The amitraz-resistance strain (RS) was screened indoors as follows: about 5000 female individuals of the SS were picked and placed on 50 flat lemon leaves, with 100 individuals per leaf. The leaves were placed into a glass dish containing sponge and water, and the four sides of the leaves were surrounded by wet cotton threads to prevent mites from escaping. The glass dish was placed under a Potter Precision Laboratory Spray Tower (Burkard Manufacturing Co., Ltd., Rickmansworth, UK), which sprayed 1 ml of amitraz emulsifiable concentrate (EC) (1:1000 dilution, 200 g l−1, Arysta LifeScience, Raleigh, NC, USA) on the mites. Then the mites were placed in an illumination incubator at 25 ± 1°C, a relative humidity of 70–80%, and a 14L:10D light cycle. Pest-free fructus citri aurantii were inoculated with the living mites 24 h after the amitraz treatment, for massive proliferation. Screening pressure was applied to establish a resistant population by spraying with amitraz EC.
Toxicity tests
The slide impregnation method recommended by the Food and Agriculture Organization of the United Nations was used for toxicity testing (Busvine, Reference Busvine1980). Based on preliminary testing, a concentration gradient of 5–7 amitraz EC doses was established, according to the equal series ratio within the mortality range of 20–90% for female adult mites. Double-sided adhesive tape was cut into 3 cm long strips and stuck to the glass slide. Healthy female adult mites were picked up with a writing brush and their back was stuck to the adhesive tape. Using this method, 30 mites were stuck to each glass slide. There were four replicates for each treatment. The slides were placed under the above-mentioned spray tower, which sprayed 1 ml of amitraz. Clear water was used for the control group. The slides were left to dry at room temperature for 15 min and then placed into an illumination incubator (temperature 25 ± 1°C, relative humidity 70–80%, 14L:10D light cycle). The dead and living mites were counted under the anatomical lens 24 h later.
Genetic analysis of the resistance
Reciprocal cross experiments
To ensure the accuracy of the experiments, female individuals of both SS and RS were collected in the quiescent stage. After molting, the unmated female mites were bred separately. RS individuals were subjected to further selection pressure indoors to remove heterozygotes. Since the female adult mites are capable of parthenogenesis to produce males, the eggs from the unmated female adult mites were collected to obtain sufficient male mites. In the direct cross experiment, 200 female SS individuals were crossed with 200 male RS individuals to obtain the F1 generation of the direct cross F1 SS♀×RS♂. In the reciprocal cross experiment, 200 female adult RS individuals were crossed with 200 male adult SS individuals to obtain the F1 generation of the reciprocal cross F1 RS♀×SS♂. The method of toxicity test introduced in section ‘Toxicity tests’ was used to calculate the dose–mortality curves for the F1 SS♀×RS♂ and F1 RS♀×SS♂ generations.
Backcross experiments
Since the female adult mites are capable of parthenogenesis to produce males, the F2 generation produced by self-cross within the F1 generation was considered the offspring by backcross between F1 generation and parents, i.e., BC(SS×RS)♀×SS♂ and BC(RS×SS)♀×RS♂. The method of toxicity test introduced in section ‘Toxicity tests’ was used to calculate the dose–mortality curves of backcross generations.
Dominant/recessive inheritance
For female mites of the F1 generation, the degree of dominance (D) was used to determine whether the inheritance of resistance trait was dominant or recessive (Stone, Reference Stone1968). The calculation formula was D = (2X 2–X 1–X 3)/(X 1–X 3), where X 1, X 2, and X 3 are the log values of LC50 for resistance in the homozygotes, heterozygotes of the F1 generation (F1 SS♀×RS♂ or F1 RS♀×SS♂), and sensitive homozygotes, respectively. The inheritance of resistance trait was considered to be completely recessive when D = −1; the inheritance of resistance trait was considered to be incompletely recessive when −1 < D < 0; the inheritance of resistance trait was considered to be incompletely dominant when 0 < D < 1; and the inheritance of resistance trait was considered to be completely dominant when D = 1.
Autosomal/cytoplasmic inheritance
If the resistance trait underwent cytoplasmic inheritance, then the trait of the offspring F1 SS♀×RS♂ and F1 RS♀×SS♂ from the reciprocal cross should be consistent with those of the female parents and show significant differences. If the resistance trait underwent autosomal inheritance, then the difference in traits should not be significant. According to David's method (David and Lehmann, Reference David and Lehmann1961), the D value of the offspring of the reciprocal cross was analyzed using analysis of variance and independent sample t-test, using the formula below:
$$\eqalign{{\rm Var}D = &\, {4 \over {{\lpar {X_1-X_3} \rpar }^2}}\bigg[{{\rm Var}\lpar {X_2} \rpar + \displaystyle{{{\lpar {X_2-X_3} \rpar }^2} \over {{\lpar {X_1-X_3} \rpar }^2}}{\rm Var}\lpar {X_1} }\rpar \cr &\quad\quad\quad\quad\quad\quad\;\,\quad\quad + \displaystyle{{{\lpar {X_2-X_1} \rpar }^2} \over {{\lpar {X_1-X_3} \rpar }^2}}{\rm Var}\lpar {X_3} \rpar \bigg]}$$
$$t = \displaystyle{{D_1-D_2} \over {\sqrt {{\rm Var}\lpar {D_1} \rpar + {\rm Var}\lpar {D_2} \rpar } }}$$The significance level was 0.05.
Monogenic/polygenic inheritance
If resistance was a monogenic trait, then the offspring BC(SS×RS)♀×SS♂ and BC(RS×SS)♀×RS♂ from the backcrosses should conform to a 1:1 separation ratio. Actual mortality, W, was plotted against expected mortality E under the same dose for the parents and backcross offspring, respectively. The difference was analyzed using a χ2 test (Preisler et al., Reference Preisler, Hoy and Robertson1990). If the difference was not statistically significant, then the two curves would coincide with each other and the inheritance was monogenic; otherwise, the two curves would not coincide with each other and the inheritance was polygenic. The formulae representing the dose–expected response E (BC) in the backcross offspring are given below:
$$E\lpar BC_{\lpar {{\rm SS} \times {\rm RS}} \rpar {\rm \female } \times {\rm SS\male }}\rpar = 0.5W_{\lpar {{\rm SS} \times {\rm RS}} \rpar {\rm \female }} + 0.5W_{{\rm SS}\male }$$
$$E\lpar {BC_{\lpar {{\rm RS} \times {\rm SS}} \rpar {\rm \female } \times {\rm RS\male }}} \rpar = 0.5W_{\lpar {{\rm RS} \times {\rm SS}} \rpar {\rm \female }} + 0.5W_{{\rm RS\male }}$$Synergism studies
To determine the metabolic capacity of three detoxification enzymes, cytochrome P450 monooxygenase, glutathione S-transferase (GST), and esterase on amitraz, the corresponding inhibitors piperomyl butoxide (PBO), diethyl maleate (DEM), and S,S,S-Tributyltrithiophosphate (DEF) were chosen and their synergistic effects with amitraz were investigated. Acetone was used to dilute PBO, DEM, and DEF to 200, 200, and 10 mg l−1, respectively. Acetone solution alone was used as the control. SS and RS individuals were treated using the slide impregnation method. Then the slides were placed into an illumination incubator. Four hours later, RS and SS individuals were treated with different concentrations of amitraz. The toxicity regression and LC50 values were calculated after 24 h of amitraz treatment. Synergism ratios (SRs) were calculated as the ratio of LC50 for amitraz without the synergist and LC50 value of amitraz with the synergists. An SR > 2 indicated a strong synergistic effect, and the higher the value, the stronger the detoxificating role of the corresponding detoxificating enzyme.
RNA extraction and cDNA library construction
RS and SS of citrus red mites (n = 200) were placed into a 1.5 ml RNase-Free centrifuge tube and immediately cryopreserved in liquid nitrogen. Each group of experiments was duplicated three times. The samples were ground in liquid nitrogen and subjected to total RNA extraction using RNA Isolator Total RNA Extraction Reagent (Vazyme, NanJing, JiangSu, China). RNA quality was determined using a Nanodrop 2000/2000C (Thermo Fisher Scientific, Waltham, MA, USA). RNA samples were purified with magnetic beads and the cDNA was synthesized with reverse transcriptase and arbitrary primers. About 200 bp DNA fragments were purified and used as templates for PCR amplification to create the cDNA library. Based on the technology of sequencing by synthesis, cDNA libraries were sequenced on Illumina HiSeq 2500.
Assembly and annotation of transcriptomes
Assembly and annotation for transcriptomes were performed according to our previous experiments (Yu et al., Reference Yu, Ding, Luo, Li, Yang, Liu, Cong and Ran2016). The reads from the sequencing platform were broken into a K-mer library by Trinity software, then these small fragments were extended from both ends and contigs were obtained. The contigs were clustered according to the paired-end information of sequences and similarity of contigs. Local assembly was conducted to generate transcripts. The longest transcript in each local region was selected for use as a unigene. All unigenes were blasted in the GO, COG, Swiss-Prot, KEGG, and NR databases for annotation.
Analysis of differential gene expression
Differential gene expression was detected using the fragments per kilobase of transcript per million mapped reads (FPKM) method (Trapnell et al., Reference Trapnell, Williams, Pertea, Mortazavi, Kwan, Van Baren, Salzberg, Wold and Pachter2010). An algorithm was implemented strictly to screen for the DEGs following the procedures described by Leng et al. (Reference Leng, Dawson, Thomson, Ruotti, Rissman, Smits, Haag, Gould, Stewart and Kendziorski2013). DEGs were defined as those with a false discovery rate (FDR) <0.01. Based on functional annotation of the DEGs, KEGG pathway and GO enrichment analyses were conducted. Pattern clustering was performed for the DEGs. The Fisher's exact test was used for the GO enrichment analysis, with a Bonferroni correction. Thus, the significantly enriched pathways and GO terms for the DEGs were obtained for further analysis.
RT-qPCR analysis
ELF1α and ACTB genes from the citrus red mites were used as internal reference genes (Niu et al., Reference Niu, Dou, Ding, Yang, Shen and Wang2012). Primers for the DEGs and internal reference genes were designed using Primer Premier 6 software (table S1). RS and SS of citrus red mites were placed into culture dishes and sprayed with 1 ml of 50 mg l−1 amitraz under a spray tower to induce gene expression. Two hours later, 200 living mites were selected for RNA extraction. Each group of experiments was duplicated three times. RNA samples from the residual susceptible strain (SS) after library construction was used as a control. HiScript II Q RT SuperMix for qPCR with gDNA wiper (Vazyme) was used for reverse transcription of extracted RNAs from the following groups: SS, SS-induced, RS, and RS-induced. ChamQ SYBR qPCR Master Mix (Vazyme) was used for real-time fluorescence quantitative PCR (RT-qPCR). The reaction was run on an ABI 7500 Fast Real-Time PCR System (Applied Biosystem, Waltham, MA, USA) and relative gene expressions were calculated using Pfaffl's method (PfaffI, Reference Pfaffl2001). The qPCR analysis was performed using 7500 software v2.0.6. Multiple comparisons were performed using one-way ANOVA (LSD method) (P < 0.05).
Result
Toxicity evaluation
Toxicity test results indicated that the LC50 of amitraz in the SS was 29.03 mg l−1 (95% CI 20.01–42.11). Six toxicity tests were performed on the female adult RS mites, which were continuously kept under selection pressure from October 2016 to March 2018 (table 1). The resistance ratio was the highest in June 2017, reaching up to 105.1, while the LC50 was 3050.97 mg l−1 (95% CI 2482.66–3749.38). In November 2017, the resistance ratio dropped to 50.24 and the LC50 was 1458.44 mg l−1 (95% CI 1188.72–1789.34). This decrease in resistance indicated fluctuating changes in the degree of resistance in the RS (fig. 1). The likely external reason was that the tolerance of both the mites and host plants was greatly challenged by the high summer temperature; the internal reason might be that the reduction in the mite population in summer led to an increased percentage of sensitive genes. Such instability of resistance trait implied that the resistant homozygotes needed to be further screened under selection pressure to ensure the accuracy of the experiment. Thus, before the genetic analysis of resistance, the optimized toxicity evaluation suggested that the resistance ratio for the RS was 81.35, while the LC50 was 2361.45 mg l−1 (95% CI 1982.89–2812.29) (March 2018).
Toxicity levels of amitraz to P. citri at different stages. The histogram indicates the indoors toxicity test result. The dotted lines represent resistance ratios compared with SS.
Amitraz toxicity evaluation in the resistance strain (RS) and the sensitive strain (SS) of P. citri
Genetic analysis of resistance
Based on reciprocal cross and backcross experiments in P. citri, toxicity tests were performed in the F1 SS♀×RS♂ and F1 RS♀×SS♂ and BC(SS×RS)♀×SS♂ and BC(RS×SS)♀×RS♂ populations, respectively (table 2). LC50 values for F1 SS♀×RS♂ and F1 RS♀×SS♂ were close (333.85 and 298.89 mg l−1, respectively). The two populations coincided within a 95% confidence interval, which indicates a comparable degree of resistance between offspring of the reciprocal crosses; the degree of resistance was intermediate between that of the sensitive parents and resistant parents. However, a great difference was observed in toxicity between BC(SS×RS)♀×SS♂ and BC(RS×SS)♀×RS♂ (195.01 and 542.94 mg l−1, respectively). The toxicity of the two backcrosses did not overlap within the 95% confidence interval, and the resistance ratio of BC(RS×SS)♀×RS♂ was about 2.78 times of that of BC(SS×RS)♀×SS♂.
Amitraz toxicity evaluation in offspring of the reciprocal cross and backcross of P. citri
n, number of mites; D, degree of dominance.
Based on the LC50 values of resistant parents, sensitive parents, and offspring of the reciprocal cross, the D value was calculated using the formula in section ‘Dominant/recessive inheritance’. For F1SS♀×RS♂ and F1 RS♀×SS♂, the D values were 0.11and 0.06, respectively (n = 630/group). Both D values were between 0 and 1, indicating that the inheritance of resistance trait was incompletely dominant. If cytoplasmic inheritance of resistance trait existed, then the trait in the F1 generation from the reciprocal cross would be influenced by the trait of the female parent and there would be a difference in the dose–mortality curve (fig. 2). Using the formula in section ‘Autosomal/cytoplasmic inheritance’, an independent sample t-test was performed for the D value of the F1 generations. If it was assumed that the variance was equal, Sig. = 0.954 > 0.05, indicating that there was no significant difference and the dose–mortality curves of the F1 generations coincided. Therefore, no cytoplasmic inheritance of the trait was involved and the resistance trait underwent chromosomal inheritance.
Dose–response curves in the F1 generation from the reciprocal cross. ‘•’ indicates the dose–response curve from sensitive parents (SS strains); ‘○’ indicates the dose–response curve from resistant parents (RS strains); ‘Δ’ indicates the dose–response curve from F1 RS♀×SS♂; ‘▾’ indicates the dose–response curve from F1 SS♀×RS♂.
Using the formula in section ‘Monogenic/polygenic inheritance’, the dose–expected response relationship was determined for the backcross offspring and a comparison was made with the actual response (fig. 3). If the resistance trait was governed by a single gene, then the difference between the expected and actual responses between the groups would not be significant. The χ2 test revealed that Σχ2 (SS×RS)♀×SS♂ = 14.44（df = 5) <χ2 0.05 = 15.21; Σχ2 (RS×SS)♀×RS♂ = 9.04（df = 6) <χ2 0.05 = 9.61. The difference between the expected and actual response of backcross offspring was below the 0.05 significance level in the χ2 test. This indicates that the resistance trait was not governed by a single gene, but probably by multiple genes. Based on the results shown above, the inheritance of amitraz resistance trait for P. citri is polygenic and incompletely dominant.
Comparison of actual and expected dose–mortality curves in the backcross offspring. ‘•’ indicates the dose–response curve from sensitive parents; ‘○’ indicates the dose–response curve from resistant parents; ‘Δ’ indicates the dose–actual response curve from BC(RS×SS)♀×RS♂; ‘▾’ indicates the dose–actual response curve from BC(SS×RS)♀×SS♂; ‘——’ indicates the dose–expected response curve from BC(RS×SS)♀×RS♂; ‘⋅⋅⋅⋅⋅⋅’ indicates the dose–expected response curve from BC(SS×RS)♀×SS♂.
Synergism studies
Toxicity tests were performed for the RS and SS of P. citri through the addition of three different synergists (table S2, table 3). After treatment with PBO, DEM, and DEF, the amitraz LC50 values for the SS were 21.37, 32.95, and 26.43 mg l−1, respectively. The SRs were 1.33, 0.86, and 1.08, respectively, all of which were below 2. Therefore, in the SS, the synergists had few synergistic effects with amitraz. This lack of synergistic effects is likely due to the high mortality induced by amitraz in SS, which makes the demonstration of synergistic effects difficult. After treatment with PBO, DEM, and DEF, the amitraz LC50 values for the RS were 711.53, 1441.51, and 906.39 mg l−1, respectively. The SRs were 2.71, 1.34, and 2.13, respectively. PBO and DEF had synergistic effects with amitraz (both had SR above 2). This indicated that after the two synergists inhibited the target detoxification enzymes, P. citri exerted a negative effect on the metabolism of amitraz. The targets of the two synergists are cytochrome P450s and esterase, respectively, which may play important roles in the metabolic detoxification of amitraz. The target of DEM is GST, which seems to have an insignificant effect on amitraz metabolism in P. citri.
Investigation on synergistic effects with PBO, DEM, and DEF in the RS and the SS of P. citri
Screening for differentially expressed genes
About 13.57 Gb of clean data were obtained through transcriptome sequencing (SRA accession number: PRJNA594464) with a Q30 ≥90.78%, indicating good sequencing quality. After assembly, 29308 unigenes were obtained from the SS and RS, with a mean length of 841.21nt and an N50 length of 1830nt. In the COG, GO, KEGG, KOG, Pfam, Swissprot, eggnog, and nr databases, 16,096 unigenes were annotated (doi:10.5061/dryad.kh1893222).
Through FPKM and Log2 (fold change) (Log2Ratio (treatment group RPKM/control group RPKM)) computation, 748 DEGs (doi:10.5061/dryad.kh1893222) were obtained from the susceptible and resistant strains of citrus red mites, including 412 upregulated genes and 336 downregulated genes (fig. 4). The gene with the highest fold of upregulation was c23816.graph_c0, with Log2FC of 8.88; the gene with the highest fold of downregulation was c11553.graph_c0, with Log2FC of −10.91. Functional annotation was performed for the DEGs. Of eight databases, 565 unigenes were annotated, including 516 from the nr database. GO terms were classified into three categories: biological process, molecular function, and cellular component. Through GO annotation, 7399 unigenes were categorized to 37,646 GO terms; among them, 240 differentially expressed unigenes were categorized as 1134 GO terms (fig. 5). Most of these GO terms participated in ‘cellular component’ (455, 40.12%), followed by ‘biological process’ (398, 35.10%) and ‘molecular function’ (281, 24.78%). The three largest sub-categories of GO terms were ‘metabolic process’ (120 GO terms), ‘catalytic activity’ (117 GO terms), and ‘binding’ (106 GO terms).
MA map of differentially expressed genes. Log2 (FPKM): logarithmic values of mean expression level; log2 (FC): logarithmic values of differentially fold change; the green dots represent downregulation of gene expression; the red dots represent upregulation of gene expression; the black dots represent genes with no significant differences in expression.
Functional annotation of differentially expressed genes based on gene ontology (GO) categorization. The abscissa represents the GO classification, the left values of the ordinate represent the percentage of the genes in the respective set, and the right values represent the number of genes.
There were 223 DEGs annotated and classified in the KEGG database, distributed over 104 pathways. The most abundant pathway was ‘lysosome’ (15 unigenes, 14.56), followed by ‘oxidative phosphorylation’ (ten unigenes, 9.71%) (fig. 6a). The pathways of ABC transporters and metabolism of xenobiotics by cytochrome P450 contain five and one unigenes, respectively; these genes may participate in the process of detoxification of acaricides in mites. Pathway enrichment was assessed using an enrichment factor. The top 20 pathways with the highest enrichment reliability are shown in fig. 6b. The pathway ‘long-term depression’ had the highest enrichment factor value, containing one gene (c11460.graph_c0) annotated serine/threonine-protein phosphatase in the nr database, with a value of Log2FC 1.23(upregulated).
KEGG classification of differentially expressed genes (a) and pathway enrichment (b). (a) The ordinate represents the different KEGG classifications and the values of abscissa mean the number of genes annotated to the pathway and their proportion to the total number of genes annotated. (b) Each dot represents a KEGG pathway, the ordinate represents the name of the pathway, and the abscissa is the enrichment factor, which indicates the ratio between the proportion of differentially expressed genes annotated to a pathway and the proportion of all genes annotated to the pathway. The bigger the enrichment factor, the more significant the enrichment level of differentially expressed genes in this pathway. The color of the circle represents the q value. The smaller the q value, the more reliable the enrichment of differentially expressed genes in this pathway. The size of the circle represents the number of genes enriched in the pathway.
The numbers of DEGs annotated in COG, KOG, and eggNOG databases were 157, 257, and 411, respectively (Fig. S1). The COG classification showed that the most DEGs were classified into ‘Carbohydrate transport and metabolism’ (30 DEGs), followed by ‘Secondary metabolites biosynthesis, transport, and catabolism’ (18 DEGs) and ‘Defense mechanisms’ (18 DEGs). The largest category in the KOG database was ‘General function prediction only’ among DEGs, followed by ‘Intracellular trafficking, secretion, and vesicular transport’ and ‘Lipid transport and metabolism’. In the eggNOG database, 185 DEGs were classified into ‘Function unknown’, only a few genes were classified into ‘Carbohydrate transport and metabolism’ (27 DEGs), ‘Signal transduction mechanisms’ (20 DEGs), and ‘Intracellular trafficking, secretion, and vesicular transport’ (21 DEGs).
Particular attention was paid to metabolism- and detoxification-related genes, such as cytochrome P450 monooxygenase (P450s), glutathione-S-transferases (GSTs), carboxyl/cholinesterases (CCEs), and ABC transporters, and 23 differentially expressed unigenes were identified in these categories (table 4). ABC transporter G family member 23-like (c23311.graph_c0) had the highest fold of upregulation, with Log2FC of 1.87; P450s 2H1-like (c24869.graph_c0) had the highest fold of downregulation, with Log2FC of −4.85. The differentially expressed unigenes included ten ABC transporters (A family, B family, F family, G family), eight P450s (4C1, 4C3, 2L1, 18a1, 6k1, etc.), four CCEs, and one dimethylaniline monooxygenase. No GSTs were found among the DEGs, which was consistent with the results of the synergism studies.
Differentially expressed unigenes likely involved in the response of P. citri to amitraz
qPCR analysis of metabolism-related genes
When stressed by amitraz, citrus red mites usually show an overexpression of metabolizing genes, which accelerate the xenobiotic metabolism and transfer toxic compounds to the extracellular domain. To verify the expression results, real-time fluorescence quantitative PCR was performed for the four SS, RS, SS-induced, and RS-induced groups. The analyses included metabolizing and detoxification enzymes and transporters. Among P450s (fig. 7), c22254.graph_c0 was annotated as P450 4C1-like in the database. After treatment with amitraz (SS-induced), the relative expression for this gene in SS was 1.03 ± 0.08 (no significant difference) compared with the control group (SS) (P < 0.05). Without amitraz treatment, the relative expression in the RS was 2.42 ± 0.13, indicating a significant upregulation compared with SS (P < 0.05); after amitraz treatment (RS-induced), the relative expression was 3.77 ± 0.30, indicating a significant upregulation compared with RS without amitraz treatment (P < 0.05). The high expression of P450 4C1-like suggested that it was involved in the detoxification of amitraz. The gene, c23144.graph_c0, was annotated as P450 2C15-like in the database. The relative expression of this gene was only 0.74 ± 0.05 in the SS-induced compared with the control group, indicating that this gene was inhibited; in the RS, the relative expression was 1.60 ± 0.05 and 1.52 ± 0.16 without and with amitraz, respectively. Thus, compared with SS, the expression of P450 2C15-like gene was upregulated, with a difference change value below 2, indicating that P450 2C15-like may play an insignificant role in amitraz metabolism. The c21051.graph_c0 and c22729.graph_c1 genes shared a similar expression pattern, indicating no differential expressions of the two genes among SS, SS-induced, and RS; however, the two genes were upregulated in RS-induced, with a change value of smaller than 2, indicating that these two P450s were insensitive to amitraz. The c21570.graph_c0 gene was annotated as dimethylaniline monooxygenase 5. In RS, this gene was inhibited. After amitraz treatment, the relative expression of this gene was upregulated to 1.71 ± 0.03, indicating a small change. Based on the results described above, of the five oxidases, P450 4C1-like was most likely to be involved in amitraz metabolism.
RT-qPCR-based verification and analysis of the differential expression of the monooxygenases, ABC transporters, and carboxyl/cholinesterases. Data are shown as means of three replicates ± standard deviation (SD); different lowercase letters indicate a significant difference between the four samples (P < 0.05); SS: susceptible strain of P. citri treated with water (control group); RS: resistant strain of P. citri treated with water; SS/RS-induced: susceptible/resistant strain of P. citri treated with amitraz.
Gene expression using qPCR was analyzed for two upregulated ABC transporters (fig. 7). In the amitraz-treated SS (SS-induced), the relative expression of c23311.graph_c0 (ABC transporters G family member 23 isoform X1) was 1.49 ± 0.07, indicating a significant upregulation compared with the control group (SS) (P < 0.05). Its relative expression in RS was 1.43 ± 0.11. After amitraz treatment (RS-induced), its relative expression was 2.55 ± 0.33. Both showed a significant upregulation, indicating that this gene may play a vital role in the transportation of metabolized toxins. The relative expressions of c22508.graph_c0 (ABC transporters G family member 23) in both SS-induced and RS were downregulated significantly (P < 0.05). There was no significant difference between the RS-induced and control group, indicating that this gene was non-vital for the transportation of amitraz metabolites.
Gene expression using qPCR was analyzed for the three upregulated CCEs (fig. 7). Compared with the control group, the expression of c17984.graph_c0 (acetylcholinesterase 4-like) did not change significantly (1.07 ± 0.06) in the SS-induced group. However, in the RS, expression of this gene was upregulated significantly (1.98 ± 0.11, P < 0.05). After amitraz treatment in RS, acetylcholinesterase 4-like was continuously upregulated (4.11 ± 1.32). Differential expressions of this gene between RS and SS indicate that this gene might be involved in the hydrolysis of amitraz. The active forms of the two esterases c20663.graph_c0 (cholinesterase 2) and c20252.graph_c0 (acetylcholinesterase 1-like) were not very active in the SS-induced, RS, and RS-induced groups, and the fold changes were small. Thus, we inferred that these two genes were insensitive to amitraz.
Discussion
If new, low-toxic acaricide products have not been developed in a timely manner, amitraz can be used as a spectral insecticide candidate to control the citrus mite for fruit farmers. However, the unreasonable application of amitraz will result in pesticide resistance in the mites. An analysis of the genetic and molecular mechanisms of amitraz resistance in mites can facilitate resistance management. According to the genetic analysis in this study, two or more genes may be responsible for amitraz resistant in citrus red mites. These genes can be categorized into metabolic-related genes and action target-related genes. Metabolic-related genes are involved in the oxidation and hydrolysis of Xenobiotic compounds via detoxification enzymes, such as cytochrome p450s, GST, and esterase. The noxious substances are transferred outside the cells via the transporters, such as ABC transporters, to reduce the toxicity of compounds on cells. Action target-related gene regulation includes the upregulation of target genes to weaken competitive inhibition or induction of target gene mutations to reduce the compound sensitivity (Van Leeuwen and Dermauw, Reference Van Leeuwen and Dermauw2016). Thus, many studies focused on metabolic-related and action target-related genes in mites once the genetic pattern of resistance was known. Van Leeuwen et al. demonstrated that the inheritance of chlorfenapyr resistance of Tetranychus urticae is polygenic incomplete recessiveness (Van Leeuwen et al., Reference Van Leeuwen, Stillatus and Tirry2004). The inheritance of spirodiclofen resistance of T. urticae is polygenic incomplete dominance (Van Pottelberge et al., Reference Van Pottelberge, Van Leeuwen, Khajehali and Tirry2009a) and the inheritance of pyridaben resistance of T. urticae is monogenic complete dominance (Van Pottelberge et al., Reference Van Pottelberge, Van Leeuwen, Nauen and Tirry2009b).
It is particularly important to know which genes are involved in amitraz metabolism and target sites in mites after understanding their genetic background. Here, synergistic experiments facilitated the screening for genes related to detoxifying enzymes in mites on the macroscopic scale. We found that cytochrome P450s and esterase of mites might take part in the process of synergistic effects on amitraz. However, the genes that play a decisive role at the molecular level were still unknown. Therefore, digital gene expression profiling analysis was applied for further screening of the candidate gene families (Snoeck et al., Reference Snoeck, Greenhalgh, Tirry, Clark, Van Leeuwen and Dermauw2017). DEGs and qPCR analysis showed that P450s, ABC transporters, and acetylcholinesterase might be involved in the detoxification of amitraz in P. citri, consistent with the results of the synergism experiments. Demaeght et al. also found that the cytochrome P450s of T. urticae were the candidate genes for spirodiclofen resistance using synergistic studies; two relevant genes, CYP392E10 and CYP392E7 were identified through DEGs analysis. Subsequently, the metabolic mechanism of CYP392E10 for spirodiclofen was confirmed through metabolic analysis in vitro and HPLC-MS (Demaeght et al., Reference Demaeght, Dermauw, Tsakireli, Khajehali, Nauen, Tirry, Vontas, Lümmen and Van Leeuwen2013).
We found ten ABC transporter genes among the DEGs of citrus red mites, including ABC-A, ABC-B, ABC-G, and ABC-F subfamilies. These ABC transporters are transmembrane proteins that participate in transporting xenobiotic metabolites out of cells, and are essential to the development of pesticide resistance (Buss and Callaghan, Reference Buss and Callaghan2008; Heckel, Reference Heckel2012). ABC transporters consist of eight subfamilies (A–H) in arthropods and 103 ABC transporters were annotated in the genome of T. urticae, most of which belong to the C, D, and H subfamilies. These proteins may be involved in acaricide detoxification in T. urticae (Van Leeuwen et al., Reference Van Leeuwen, Vontas, Tsagkarakou, Dermauw and Tirry2010; Dermauw et al., Reference Dermauw, Osborne, Clark, Grbić, Tirry and Van Leeuwen2013). In addition, mutations in the SeABCC2 gene may be responsible for the resistance of Spodoptera exigua to Bt toxin (Park et al., Reference Park, González-Martínez, Navarro-Cerrillo, Chakroun, Kim, Ziarsolo, Blanca, Cañizares, Ferré and Herrero2014). Among the DEGs in P. citri, the PcABCG gene was upregulated in the resistant strain of mite. The PcABCG gene may be involved in amitraz detoxification in mites. Current reports indicate that the G subfamily of transporters is involved in insect pigmentation, molting, and reproduction; however, involvement in acaricide resistance has not been demonstrated (Broehan et al., Reference Broehan, Kroeger, Lorenzen and Merzendorfer2013).
Conventional studies on resistance mechanisms are mainly based on the investigation of synergistic effects, enzymatic activity determination, or real-time fluorescence quantitative PCR to analyze the metabolic resistance. Suspected target gene mutations are subjected to base alignment through Sanger sequencing in the resistant and sensitive strains. However, these methods have the defect of being researcher-dominated and it is difficult to objectively identify real resistance genes at the whole genome and transcriptome level. In the era of the genome, bulked segregant analysis (BSA) and quantitative trait locus (QTL) resequencing technologies can overcome this investigator bias. BSA is favored for monogenic inheritance, while QTL and GWAS are recommended for polygenic inheritance. BSA also applies to polygenic inheritance, such as host adaptation studies (Van Leeuwen and Dermauw, Reference Van Leeuwen and Dermauw2016). By performing genetic analysis of resistance, Van Leeuwen et al. found that etoxazole resistance of T. urticae was a monogenic recessive trait. By crossing the extremely sensitive and extremely resistant strains, the trait inheritance population and the resistant strain were screened from the F6 generation for using a bulked segregant mapping analysis. The I1017F sense mutation in the chitin synthase 1 gene was identified by this method (Van Leeuwen et al., Reference Van Leeuwen, Demaeght, Osborne, Dermauw, Gohlke, Nauen, Grbić, Tirry, Merzendorfer and Clark2012). Demaeght et al. also used bulked segregant analysis to confirm that in T. urticae, the mutations related to resistance to the mite growth inhibitors, clofentezine, hexythiazox, and etoxazole, were all I1017F in the chitin synthase 1 gene (Demaeght et al., Reference Demaeght, Osborne, Odman-Naresh, Grbić, Nauen, Merzendorfer, Clark and Van Leeuwen2014). Jagadeesan et al. established the crosses between phosphine-sensitive and phosphine-resistant Tribolium castaneum populations. Then low-coverage sequencing was applied to the DNA pool from the offspring with trait separation. Candidate regions and relevant SNPs were identified, hence narrowing the scope of target genes (Jagadeesan et al., Reference Jagadeesan, Fotheringham, Ebert and Schlipalius2013). Park et al. found that ABCC transporters were involved in the resistance of S. exigua to Bt toxins through bulk segregant analysis (Park et al., Reference Park, González-Martínez, Navarro-Cerrillo, Chakroun, Kim, Ziarsolo, Blanca, Cañizares, Ferré and Herrero2014). To study non-monogenic resistant trait, Snoeck et al. applied high-resolution QTL mapping to localize target genes related to the resistance of T. urticae to METI-I acaricides; a PSST homologue, cytochrome P450 reductase, and nuclear hormone receptor 96 were identified (Snoeck et al., Reference Snoeck, Kurlovs, Bajda, Feyereisen, Greenhalgh, Villacis-Perez, Kosterlitz, Dermauw, Clark and Van Leeuwen2019). In the future, whole genome resequencing for P. citri will be performed to localize resistant genes response to amitraz in mites.
Conclusion
The genetic analyses for the sensitive and the amitraz-resistant strains of P. citri indicate that the inheritance of resistance trait in mite is incomplete dominant, governed by polygenes on the chromosome. Synergism studies demonstrated that cytochrome P450s and esterase may play important roles in the detoxification of amitraz. Based on differential gene analysis, 23 metabolism-related genes of P. citri were identified, consistent with the results of the synergism studies. Real-time PCR verification implied that P450s, ABC transporters, and acetylcholinesterase might influence the detoxification of amitraz by P. citri. These results provide the genetic and molecular foundation for the management of pest mite resistance.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485320000267
Financial support
This research was supported by the National Key Research and Development Program (2018YFD0201508); the National Agricultural Innovation Project (CAAS-XTCX-2016013); the Chongqing Scientific Research Project (cstc2016shms-ztzx80003, cstc2017shms-xdny0443); the Fundamental Research Funds for the Central Universities (XDJK2015C014); and the Research Fund for the Doctoral Program of Southwest University (SWU116061).
