Comparative Transcriptome Analysis Reveals the Mechanism Related to Fluazinam Stress of Panonychus citri (Acarina: Tetranychidae)

Simple Summary The citrus red mite, Panonychus citri, is an important pest that causes serious citrus production losses in China. The insecticide fluazinam has a good control effect on the pest mites; however, its mechanism of action on mites remains unclear. In this study, we analyzed the transcriptomic sequencing and differential expression genes in P. citri treated with fluazinam, and identified some of the genes potential involved in detoxification metabolism related with the fluazinam exposure. Evaluating the efficacy of fluazinam, and analyzing the transcriptome data of P. citri under fluazinam stress, potentially provide a new agent for prevention and control of P. citri, and also preliminary research results for exploring the mechanism of action of fluazinam on P. citri. Given the up-regulated expression levels of genes for Mn-superoxide dismutase and catalase, we speculate that they play an important role in fluazinam-stress action on P. citri. Abstract The use of a large number of chemical acaricides to control these pest mites has led to an increasing problem of pesticide resistance, which has always been the difficulty in integrated pest management (IPM). Fluazinam has a good control effect on Panonychus citri, the serious pest on citrus; however, we only know the mechanism of action of fluazinam as a fungicide and its mechanism of action on mites remains unclear. Through analysis using Illumina high-throughput transcriptomic sequencing and differential expression genes in P. citri treated with fluazinam, 59 cytochrome P450 genes, 23 glutathione s-transferase genes, five carboxylate esterase genes, 11 superoxide dismutase genes and 15 catalase genes were identified. The Gene Ontology enrichment and the enrichment of KEGG results showed that the treatment were enrichment for redox enzyme pathways. Evaluating the efficacy of fluazinam, and analyzing the transcriptome data of P. citri under fluazinam stress, potentially provide a new agent for prevention and control of P. citri, and also preliminary research results for exploring the mechanism of action of fluazinam on P. citri. Given the up-regulated expression levels of genes for Mn-superoxide dismutase and catalase, we speculate that they play an important role in fluazinam-stress action on P. citri.


Introduction
Panonychus citri (Acarina: Tetranychidae), also known as citrus red mite, is an important fruit tree pest in many citrus-producing areas of China, and can harm many kinds of citrus plants. Chemical control is still the main method to control citrus red mite because of convenience,

P. citri Samples
The P. citri were collected from a citrus orchard in Zhejiang Province, and raised on citrus seedlings in a greenhouse in Zhejiang A&F University (25 ± 1 • C, 80 ± 5% RH, 14 h:10 h light:dark). Fluazinam TC was dissolved in acetone solution, and then distilled water was used to prepare the solution of sub-lethal concentration. Samples of active P. citri were collected at 6, 24, and 48 h after treatments, with about 300 individuals in each sample. The RNA extraction was then performed for subsequent transcriptome sequencing [25]. Acetone treatment at the same concentration was used as the control. All treatments were replicated three times. All RNA samples were immediately stored in a refrigerator at −80 • C.

RNA Extraction and Library Construction
The total RNA extraction kit of Trizol (UNlQ-10 Column Trizol Total RNA Isolation Kit, Sangon Biotech, Shanghai, China) was used. Sequencing libraries were generated using NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (New England BioLabs, Ipswich, MA, USA) according to manufacturer's protocol. Briefly, magnetic beads with Oligo(dT) were used to concentrate the mRNA. Then fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer. First strand of cDNA was synthesized with random hexamers and M-MuLV Reverse Transcriptase. Both DNA polymerase I and RNase H were used to synthesize the double-stranded cDNA, and then the double-stranded cDNA was purified by AMPure XP beads. The end of the obtained double-stranded cDNA was repaired, a tail was added and the sequencing joint was connected. Then AMPure XP beads were used to screen the size of the segments (preferentially 250-300 bp). Finally, PCR amplification and purification were performed to obtain the final library. Qubit ® 2.0 was used for initial quantification and dilution of the library, followed by determination of the insert size of the library. Illumina HiSeq sequencing was performed after the library inspection and 150 bp paired-end reads were generated by Novogene (Beijing, China). The transcriptome raw data has been released already with ID: PRJNA669340.

Differential Expression Analysis
The input data of differential gene expression is the read-count data obtained in analysis of gene expression level. We used DESeq R package (1.10.1) for analysis [28], which provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution, and the screening threshold was Padj < 0.05 and |log 2 FoldChange| > 1. In the process of software differences analysis, we corrected the p-value using Benjamini and Hochberg's approach obtained from the original hypothesis test.

Enrichment and KEGG Enrichment Analysis
The GO enrichment analysis method was using GOseq [29] R package, which is based on Wallenius non-central hypergeometric distribution, and can adjust for gene length bias in DEGs. The KEGG database is a resource for understanding high-level functions and utilities of biological systems from molecular-level information [30] (http://www.genome.jp/kegg/). Significant pathway enrichment analysis was conducted with KEGG pathways as units, and hypergeometric tests were used to identify pathways with significant enrichment of differentially expressed genes relative to all annotated genes by using KOBAS 2.0 software [31].

Splicing Transcript Length Distribution
Transcriptional sequences obtained from Trinity were used as reference sequences for subsequent analysis. The longest cluster sequence was obtained by Corset hierarchical clustering for subsequent analysis. The lengths of transcripts and clustering sequences were statistically analyzed (Table 1). After splicing, the total number of genes was 28,831, of which 10,492 were 200-500 bp in length, accounting for 36.39% of the total number of genes. There were 8395 genes with a length of 500-1000 bp, accounting for 29.12% of the total number; 4151 genes of length 1000-2000 bp accounted for 14.39%, and 5793 genes of length exceeding 2000 bp accounted for 20.10%.

Gene Annotation Success Rate Statistics
In order to provide comprehensive gene function information, the obtained unigenes were annotated with gene function in seven databases, and the success rate of annotation in the seven databases was calculated ( Table 2). In total, 18,771 unigenes (65.1%) were annotated in at least one database.

Unigene Homology Analysis
Homology analysis of unigenes successfully annotated in NR showed that homologs of these genes included multiple species ( Figure S1). Of unigenes in the citrus red mite, 48.7% had the highest homology with genes related to the two-spotted spider mite (Tetranychus urticae Koch), followed by Insects 2020, 11, 730 5 of 14 the brown rice planthopper (Nilaparvata lugens (Stal)) with 3.3%. In addition, Hydra (Anthoathecatae, Hydridae) and western predatory mite (Galendromus occidentalis) had homology with about 2.0%, and 1.6% showed high homology with an amoeba, Capsaspora owczarzaki.

GO Function Analysis
After GO annotation of genes, the genes that were annotated successfully were classified according to the next level of the three major categories of GO ( Figure S2), Among them, unigenes were divided into cell components, followed by biological processes and finally molecular functions. The three functional groups with the most unigene annotations were, in turn, 8591 unigenes in the cellular processes of the cellular components. There were 8061 unigenes in metabolic process and 7276 unigenes in the bindings in the molecular function. Note that the three functional groups with the least unigenes follow: One unigene in nucleoid in biological process, one unigene in cell aggregation in cell component, and two unigenes in extracellular matrix component in biological process.

KOG Function Analysis
KOG, euKaryotic Ortholog Groups, is a database of eukaryotic organisms based on lineal homologous relationship of genes. The obtained P. citri unigenes were compared to the KOG database. There were 8601 successfully annotated unigenes ( Figure S3). These unigenes were assigned to 25 different functional families in the KOG database. The three functional families with the most annotated unigenes were "Translation, ribosomal structure, and biogenesis", and "General function prediction only" and "Posttranslational modification, protein turnover, chaperones (post-translation modification, protein conversion, chaperones)". These families contained 1247, 1116, and 1048 unigenes, respectively, accounting for 14.49, 12.98, and 12.18% of successfully annotated genes. Note the functional families with the least unigenes were "Defense mechanisms" with 44, "Nuclear structure" with 26 and "Cell motility" with eight, accounting for 0.51, 0.3, and 0.09% of successfully annotated genes, respectively ( Figure S3).

Differential Gene Analysis
Based on the statistics of differentially expressed genes between different treatments and control, a quantitative relationship between the number of differential genes and the total number of genes was obtained. After 6 h of treatment, only 28 differential genes were up-regulated and 487 were down-regulated ( Figure 1); after 24 h, 368 were up-regulated and 1997 down-regulated and, after 48 h, 191 were up-regulated and 1381 down-regulated. After 24 and 48 h of treatment, the number of differentially expressed genes was significantly higher than that after 6 h, and the proportion relative to the total number of genes increased significantly. The number of differentially expressed genes was the highest at 24 h after treatment. The number of down-regulated genes in the total gene composition of the three treatment groups, significantly more genes were down-regulated than up-regulated (differential gene screening conditions were |Log2 FoldChange| > 1 and Padj < 0.05).
h, 191 were up-regulated and 1381 down-regulated. After 24 and 48 h of treatment, the number of differentially expressed genes was significantly higher than that after 6 h, and the proportion relative to the total number of genes increased significantly. The number of differentially expressed genes was the highest at 24 h after treatment. The number of down-regulated genes in the total gene composition of the three treatment groups, significantly more genes were down-regulated than upregulated (differential gene screening conditions were |Log2 FoldChange| > 1 and Padj < 0.05).

GO Enriched Gene Analysis
The GO gene enrichment analysis is the application of biological information databases and statistical tools to assign enriched target genes to biological pathways or modules of known functions, thereby allowing in-depth study of gene function from a biological perspective. The purpose is to screen out two groups or enriched gene sets with different expression levels between multiple groups. The 10 functional options with the most significant GO enrichment in different treatment groups were selected for statistical comparison and the results are shown in Table 3. At 6 h, GO enrichment contained the largest number of differential genes: 255 of metabolic process (GO: 0008152) followed by 209 of catalytic activity (GO: 0003824). The most enriched function was fungal-type cell wall (GO: 0009277); however, only six of the differential genes were related to this option. In addition, oxidation-reduction process (GO: 0055114) and oxidoreductase activity (GO: 0016491) also showed some degree of enrichment.
Based on the enrichment of GO function at this time point, we speculate that within a short period of application of fluazinam, the stress response of P. citri to fluazinam was manifest by significant changes in metabolic activity and catalytic activity in the body. At the same time, the body's redox process and energy metabolism were also somewhat affected. At 24 h, GO enrichment statistics showed that the functional item with the highest enrichment level was oxidation-reduction process, followed by oxidation activity, which included the functional item (metabolic process) with the largest number of differential genes. A total of 1059 differential genes were related to the function of metabolic process. Based on the 24-h GO enrichment, we speculate that significant changes in redox process and energy metabolism were closely related to the mechanism of action of fluazinam on P. citri. In addition, the abnormal metabolic function may suggest the subsequent response of P. citri to fluazinam. The functional item with the highest enrichment level at 48 h was cofactor binding (GO: 0048037), followed by coenzyme binding (GO: 0005319) and oxidation activity (GO: 0016491). Among the 48-h statistical GO enrichment functions, redox was still an important item of enrichment, and cofactor coenzyme binding also showed a high degree of enrichment.

KEGG Enrichment Pathway Analysis
The results of the KEGG enrichment pathway analysis are shown in Figure 2. Compared with GO enrichment results, at 6 h, the KEGG enrichment pathway information showed that the significant results were degradation of valine, leucine and isoleucine degradation, as well as on the ribosome. At 24 h, the KEGG enrichment pathways were proteasome and peroxisome, similar to the abnormal results of redox and energy metabolism found in the GO enrichment function at 24 h. At 48 h, the KEGG enrichment result showed that the enrichment pathway was glycolysis/gluconeogenesis and citrate cycle (TCA cycle) still showed a close relationship with energy synthesis and metabolism in the living body. The GO and KEGG enrichment results indicate that the mechanism of action of fluazinam on P. citri was important and associated with its energy metabolism.

KEGG Enrichment Pathway Analysis
The results of the KEGG enrichment pathway analysis are shown in Figure 2. Compared with GO enrichment results, at 6 h, the KEGG enrichment pathway information showed that the significant results were degradation of valine, leucine and isoleucine degradation, as well as on the ribosome. At 24 h, the KEGG enrichment pathways were proteasome and peroxisome, similar to the abnormal results of redox and energy metabolism found in the GO enrichment function at 24 h. At 48 h, the KEGG enrichment result showed that the enrichment pathway was glycolysis/gluconeogenesis and citrate cycle (TCA cycle) still showed a close relationship with energy synthesis and metabolism in the living body. The GO and KEGG enrichment results indicate that the mechanism of action of fluazinam on P. citri was important and associated with its energy metabolism.

Genetic Verification
In order to verify the authenticity of the transcriptome data, a total of 10 genes were randomly selected from three different treatments (6, 24 and 48 h) for qPCR experiments. The overall gene expression trend was consistent with the transcriptome data except two genes (Table 4).

Genetic Verification
In order to verify the authenticity of the transcriptome data, a total of 10 genes were randomly selected from three different treatments (6, 24 and 48 h) for qPCR experiments. The overall gene expression trend was consistent with the transcriptome data except two genes (Table 4).

Discussion
Cytochrome P450, glutathione s-transferase and carboxylate esterase are three important supergene families in the metabolic resistance of insects (and pest mites) [32][33][34]. Ding et al. indicated that P450 gene expression was over-transcribed in acaricide exposures of P. citri. The gene CYP4CF1 is induced by pyridaben and CYP4CL2 expression is induced by abamectin, azocyclotin, and pyridaben [35]. Transcriptome and differential expression analysis of Tetranychus cinnabarinus following treatment with a sub-lethal concentration of β-sitosterol revealed that the detoxification genes of Tetranychus cinnabarinus were up-regulated, including 28 carboxyl/cholinesterases and ABC transporter C-class. Some defense-related proteins were also activated, such as toll-like receptors, legumain and serine proteases [36]. In this study, a total of 59 P450-related genes were correlated with fluazinam, among which the CYP2, CYP3, and CYP4 clans were the most widely distributed. These genes are involved in the metabolic pathways of P. citri to exogenous substances (ko01100). In addition, 23 glutathione-S-transferase genes and five carboxylesterase genes were identified among the differentially expressed genes. Pathway annotations suggested that these genes might be involved in the metabolic process of P. citri in response to exogenous chemicals (ko00982 and ko00983). Therefore, we hypothesize that these genes are closely related to fluazinam stress of P. citri.
Many researchers have carried out "omics" studies on mites, including reproduction and sexual differentiation, as well as the detoxification and resistance mechanisms related to acaricides [37,38]. Previous authors considered that the analysis of pest mites at the molecular level based on high-throughput sequencing was helpful to identify the genes related to xenobiotic metabolism [36] and insecticide resistance [39]. Through transcriptome data analysis of Tetranychus cinnabarinus, 10 gene categories related to pesticide resistance were identified: P450, GST, CarE, AChE, GluCl, nAChR, GABA receptor, sodium channel, ATPase, and Cyt b genes [39]. Further GO enrichment analysis using data on differentially expressed genes of fenpropathrin-susceptible and -resistant strains of Tetranychus cinnabarinus showed that the differences in catalytic activity and binding between the two strains were the most significant, suggesting that these two biological functions might affect the sensitivity of Tetranychus cinnabarinus to fenpropathrin [40]. Comparative transcriptome study of resistant and susceptible strains of P. citri resulted in identification of 211 metabolism genes and target genes related to insecticide resistance [22]. These conclusions provided a basis and reference for our research. In our study, catalytic activity was significantly enriched and higher than for other GO functional items at 6 h after fluazinam treatment. Thus, it is likely that catalytic activity of mites in vivo changed rapidly after being stimulated by exogenous chemical fluazinam. The genes enriched in this functional term may also be involved in fluazinam-stress in P. citri.
In addition, GO and KEGG enrichment analyses showed that the influence of fluazinam on P. citri was not only in catalytic activity and metabolic pathway, but also in energy synthesis and redox process. The mechanism of action of fluazinam in other organisms is uncoupling and blocking energy synthesis [41,42]. In this study, a total of 11 superoxide dismutase (SOD) and 15 catalase genes were identified. Our results revealed SOD and catalase genes involved in the mechanism of action of P. citri on exogenous fluazinam. The SOD is an active enzyme existing widely in animals, plants and microorganisms, which can remove the harmful substances produced in biological metabolism in the innate immune system [43,44]. The SOD can catalyze the disproportionation of superoxide free-radicals in organisms and eliminate these free radicals. Moreover, SOD can block the damage to cells caused by oxygen free-radicals and repair the damaged cells in time [45]. The SODs can be divided into three types: Cu-ZnSOD, MnSOD, and FeSOD [46]. The Cu-ZnSODs mainly exist in the cytoplasm of eukaryotic cells and are one of the most widely distributed SODs in nature; MnSODs are mainly concentrated in the mitochondria and prokaryotic cells of eukaryotic cells, and FeSODs mainly exist in prokaryotic cells [44]. Much recent research has shown that SOD plays a major role in protecting insects against various environmental stresses [47], including insecticides. Previous studies showed that SOD plays a major role in the process of insect (and mite) response to various adverse environments, including pesticide [48] and temperature stresses [49,50].
The gene encoding a MnSOD (PcSOD3) plays a role in the anti-oxidative damage to P. citri. Gene cloning and functional exploration of the SOD of P. citri showed that MnSOD (PcSOD3) in P. citri has more strong antioxidant ability than Cu-ZnSOD (PcSOD1 and PcSOD2) [51]. In addition, they also found that after P. citri was subjected to external environmental stresses such as acaricides, extreme temperature and ultraviolet radiation, three different types of SOD genes in the body showed up-regulation to varying degrees, with the gene for MnSOD (PcSOD3) greatly up-regulated, suggesting that it plays an important role in P. citri response to adverse environmental stress [51,52]. Treatment of zebrafish with a half lethal concentration of fluazinam resulted in no significant change in expression levels of genes for Cu-ZnSOD and catalase (CAT), while that for MnSOD significantly increased [53]. In this study, genes for MnSOD and CAT in P. citri were significantly up-regulated by fluazinam treatment. The enzyme MnSOD mainly exists in mitochondria and plays an antioxidant role, whereas Cu-ZnSOD exists outside mitochondria. Based on the experimental data, we speculate that fluazinam may promote the generation of a large number of oxidative free-radicals in the mitochondria of P. citri, leading to abnormal mitochondrial function and thus playing a role in killing the mite. Oxidative stress is a key factor in the mode of action of pesticide imidacloprid at low doses to Drosophila. This pesticide induced increase levels of reactive oxygen species (ROS) to affect mitochondrial function, energy levels, the lipid environment, and transcriptomic profiles [54].

Conclusions
In this study, 59 cytochrome P450 genes, 23 glutathione s-transferase genes, five carboxylate esterase genes, 11 superoxide dismutase genes and 15 catalase genes related with the fluazinam exposure were identified by transcriptomic analysis and differential expression genes in P. citri treated with fluazinam. Given the up-regulated expression levels of genes for Mn-superoxide dismutase and catalase, we speculate that they play an important role in the mechanism of fluazinam action on P. citri. Based on the experimental data, we speculate that fluazinam may promote the generation of a large number of oxidative free-radicals in the mitochondria of P. citri, leading to abnormal mitochondrial function and thus playing a role in killing the mite, that they play an important role in the mechanism of fluazinam action on P. citri. The transcriptome data of P. citri are important in revealing the mechanism of acaricidal action of fluazinam, and provide theoretical basis for the control of resistance of mite in IPM.