Transcriptome Analysis and Identification of Insecticide Tolerance-Related Genes after Exposure to Insecticide in Sitobion avenae

Aphids cause serious losses to the production of wheat. The grain aphid, Sitobion avenae, which is the dominant species of aphid in all wheat regions of China, is resistant to a variety of insecticides, including imidacloprid and chlorpyrifos. However, the resistance and mechanism of insecticide tolerance of S. avenae are still unclear. Therefore, this study employed transcriptome analysis to compare the expression patterns of stress response genes under imidacloprid and chlorpyrifos treatment for 15 min, 3 h, and 36 h of exposure. S. avenae adult transcriptome was assembled and characterized first, after which samples treated with insecticides for different lengths of time were compared with control samples, which revealed 60–2267 differentially expressed unigenes (DEUs). Among these DEUs, 31–790 unigenes were classified into 66–786 categories of gene ontology (GO) functional groups, and 24–760 DEUs could be mapped into 54–268 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Finally, 11 insecticide-tolerance-related unigenes were chosen to confirm the relative expression by quantitative real-time polymerase chain reaction (qRT-PCR) in each treatment. Most of the results between qRT-PCR and RNA sequencing (RNA-Seq) are well-established. The results presented herein will facilitate molecular research investigating insecticide resistance in S. avenae, as well as in other wheat aphids.


Introduction
Wheat, Triticum aestivum, is considered one of the most important cereals in China, as well as around the world. Aphids cause losses of more than 10% of the harvest in an average year and over 30% in years with serious damage [1]. The grain aphid, Sitobion avenae, is the dominant aphid in all wheat regions of China and is also responsible for the greatest loss to wheat production [2]. Both the aphid adult and nymph feed on the wheat leaves, stems, and ears, causing the plants to stunt, be unable to ear, or even die. Aphids also damage wheat by the transmission of plant viruses, such as barley yellow dwarf virus, and impact photosynthesis by the production of honeydew, all of which lead to wheat yield loses and poor quality [1][2][3].
In China, insecticides treatment is the most prevalent management strategy used against grain aphids when planting wheat. Identifying the specific genes involved in insecticide detoxification and

Sequencing Data Analysis
SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) with default parameters were used to trim and control the quality of raw paired reads. The clean data of insecticide-treated and control samples were used to conduct de novo assembly with Trinity (http://trinityrnaseq.sourceforge.net/) [12]. The assembled transcripts were searched against the NCBI protein nonredundant (NR), Swiss Prot, Pfam, and Cluster of orthologous groups (COG) databases using BLASTX, the proteins that had the highest sequence similarity with the given transcripts were used to retrieve their function annotations, and typical cut off E values less than 1.0×10 −5 were set. Gene ontology (GO) annotations of unique assembled transcripts for describing biological processes, molecular functions, and cellular components were generated using BLAST2GO (https://www.blast2go.com/) [13]. The Kyoto Encyclopedia of Genes and Genomes (KEGG, http: //www.genome.jp/kegg/) [14] was used to perform metabolic pathway analysis. The raw reads were deposited in the NCBI Sequence Read Archive database under accession numbers SRR8953735-8953755.

Analysis of Differential Expressed Unigenes (DEUs), GO Annotation, and Pathway Enrichment
Gene expressed levels were assessed by RSEM (RNA-Seq by Expectation Maximization) [15] for each sample using the basic method of clean data mapping back to the assembled transcriptome to obtain the read count of each gene. The DESeq2 was used to perform the DEUs' analysis of control groups and insecticide treatment groups [16]. The relative expression of unigenes was calculated by dividing the unigene's fragments per kilobase per million mapped reads (FPKM) value into insecticide-treated groups with the same unigene FPKM value in control groups. The R statistical package software EdgeR (Empirical analysis of Digital Gene Expression in R, http://www.bioconductor. org/packages/2.12/bioc/html/edgeR.html) was utilized for differential expressed analysis, and the R code for detect differential expression of unigenes between treated with imidacloprid for 15 min and control was provided in the Supplementary File S1 [17]. In addition, functional enrichment analysis, including GO and KEGG, was performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways at Bonferroni corrected p-values < 0.05 when compared with the total transcriptome background. GO functional enrichment and KEGG pathway analysis were conducted using Goatools (https://github.com/tanghaibao/Goatools) and KOBAS (http://kobas.cbi.pku.edu.cn/) [18].

Sequence Confirmation and qRT-PCR Validation
Eleven DEUs assembled sequences, including one ABC transporter, one glutathione s-transferases (GST), one esterase, three cytochrome P450, two uridine 5 -diphospho-glucuronosyltransferase (UGT), three trypsin, and one reference gene nicotinamide adenine dinucleotide (NADH) were confirmed and the irrelative expressions validated by qRT-PCR. Primer premier 5.0 was used to design specific primers to confirm the assembled sequences, and then reverse transcription PCR was conducted (Supplementary File S2). PCR products were then analyzed by gel electrophoresis, and expected DNA bands were extracted using an Agarose Gel Extraction Kit (Takara Biotechnology (Dalian) Co., Ltd., Dalian, Liaoning, China). Then DNA was sub-cloned into the vector (pEASY-T1 Simple Cloning Kit, Beijing TransGen Biotech Co., Ltd., Beijing, China) according to the manufacturer's protocols and a 3730 DNA analyzer was used to determine the nucleotide sequences. The confirmed sequences were deposited in the GenBank database under accession numbers MN481369 to MN481380.The relative expression of 11 confirmed sequences was validated by qRT-PCR, while NADH was chosen as the reference gene [19]. Beacon Designer 7.0 was used to design the specific primers for qRT-PCR (Supplementary File S2). One microgram total RNA was used to synthesize the cDNA after removal of genomic DNA (PrimeScript™ RT Reagent Kit with gDNA Eraser, Takara Biotechnology (Dalian) Co., Ltd., Dalian, Liaoning, China). SYBR Premix (Takara Biotechnology (Dalian) Co., Ltd., Dalian, Liaoning, China) was used for qRT-PCR conduct on a CFX96 Real-Time PCR Detection System (BioRad Laboratories, Inc., Hercules, CA, USA). The thermal cycling conditions were 45 cycles of 95 • C for 30 s, 58 • C for 30 s, and 72 • C for 30 s. The qRT-PCR was done accordingly to MIQE (minimum information for Q-PCR experiment) [20]. The amplification efficiency of each pair of specific primers was calculated (Supplementary File S3) [21]. Three biological samples, with two technical replicates, were used to determine the relative expression of unigenes. The NCBI accession numbers of confirmed sequences and the primers designed for confirmation and qRT-PCR are listed in File S2.

Bioassay, RNA Sequencing, Assembly, and Annotation
The bioassay results are shown in Table 1  ) unigenes were annotated in the NR, Swiss-Prot, Pfam, COG, GO, and KEGG databases, respectively. The top three functions among the GO annotations were cellular processes (15.30%), metabolic processes (13.82%), and single organism processes (11.76%) ( Figure 1A). In the KEGG mapped results, the pathways were divided into six main categories, of which human diseases (10,915), organizational (8115), and metabolism (8067) accounted for the top three ( Figure 1B, Supplementary File S7). In the COG, the largest number of annotations (903 unigenes) were translation, ribosomal structure, and biogenesis, which belonged to information storage and processing, followed by 521 unigenes annotated to posttranslational modification, protein turnover, and chaperones, which belonged to cellular process and signaling, and the third was 323 unigenes annotated to general function prediction only, which has been poorly characterized ( Figure 1C, Supplementary File S8).

Unigene Expression Analysis for Insecticide Treatments
The identification of DEUs of insecticide-treated samples was compared with the control based on the FPKM value. Unigenes were considered DEUs only if the fold change (FC) expression ratios of insecticide-treated samples versus the control sample were larger than two or less than 0.5. A total of 60-2267 unigenes were considered DEUs after pesticide treatment ( Figure 2). The number of DEUs increased as the insecticide treatment time increased. S. avenae treated with chlorpyrifos for 36 h had the most downregulated unigenes (1213), while the minimum number of downregulated unigenes (47) Figure  3B). However, for the same time duration, different insecticides treatments also had intersection DEUs. After 15 min, different insecticides treatment had 20 intersection DEUs ( Figure 3C), while for 3 h, there were 103 intersection DEUs when treated with imidacloprid or chlorpyrifos ( Figure 3D), and for 36 h, there were 1518 intersection DEUs when treated with imidacloprid or chlorpyrifos ( Figure 3E).

GO Classification and KEGG Pathway Identification of DEUs
DEUs generated from different comparisons were assigned into GO term analysis. A total of 31-790 unigenes were classified into 66-786 categories of GO functional groups ( Figure 4). The DEUs of imidacloprid treated for 36 h in comparison with the control had the largest number of unigenes that could be classified into GO terms, of which 790 unigenes were assigned into 786 GO terms. The DEUs of imidacloprid treated for 15 min had the lowest number of unigenes that could be assigned into GO terms, of which only 31 were categorized into 66 GO terms (Supplementary File S10). To identify possible active biological pathways of DEUs, the unigenes were mapped into the KEGG pathways. In the DEUs of different comparisons, 24−760 DEUs could be mapped into 54−268 KEGG pathways ( Figure 5). The imidacloprid treatment for 36 h had the largest DEUs (760) mapped into the largest

GO Classification and KEGG Pathway Identification of DEUs
DEUs generated from different comparisons were assigned into GO term analysis. A total of 31-790 unigenes were classified into 66-786 categories of GO functional groups ( Figure 4). The DEUs of imidacloprid treated for 36 h in comparison with the control had the largest number of unigenes that could be classified into GO terms, of which 790 unigenes were assigned into 786 GO terms. The DEUs of imidacloprid treated for 15 min had the lowest number of unigenes that could be assigned into GO terms, of which only 31 were categorized into 66 GO terms (Supplementary File S10). To identify possible active biological pathways of DEUs, the unigenes were mapped into the KEGG pathways. In the DEUs of different comparisons, 24−760 DEUs could be mapped into 54−268 KEGG pathways ( Figure 5). The imidacloprid treatment for 36 h had the largest DEUs (760) mapped into the largest number of pathways (268), while imidacloprid treatment for 15 min led to the lowest number of DEUs (24) mapped into the least pathways (54) (Supplementary File S11). The ten most up-and downregulated DEUs after imidacloprid and chlorpyrifos treatment for 36 h compared with the control are listed in Tables 2 and 3

Insecticide Tolerance Related Unigenes Analysis and Quantitative Real-Time PCR (qRT-PCR) Validation
The DEUs of insecticide-treated and control samples that were considered to be insecticide tolerance-related unigenes are listed in Table 4. As the insecticide treatment time increased, the number of insecticide tolerance-related DEUs increased. After treatment with insecticides for 36 h, cuticle protein had the largest number of DEUs, followed by the ABC transporter and trypsin. After insecticide treatment for 36 h, the same 39 DEUs related to cuticle proteins were found in both the imidacloprid and chlorpyrifos treatment, while there were 50 DEUs in the imidacloprid treatment

Insecticide Tolerance Related Unigenes Analysis and Quantitative Real-Time PCR (qRT-PCR) Validation
The DEUs of insecticide-treated and control samples that were considered to be insecticide tolerance-related unigenes are listed in Table 4. As the insecticide treatment time increased, the number of insecticide tolerance-related DEUs increased. After treatment with insecticides for 36 h, cuticle protein had the largest number of DEUs, followed by the ABC transporter and trypsin. After insecticide treatment for 36 h, the same 39 DEUs related to cuticle proteins were found in both the imidacloprid and chlorpyrifos treatment, while there were 50 DEUs in the imidacloprid treatment (four upregulated and 46 downregulated) and 43 DEUs in the chlorpyrifos treatment (three upregulated and 40 downregulated). ABC transporter possessed the second largest number of DEUs after insecticide treatment for 36 h (29 DEUs), with eight upregulated and 21 downregulated in the imidacloprid treatment and 25 DEUs with eight upregulated and 17 downregulated in the chlorpyrifos treatment. Among these DEUs, the same 19 DEUs related to ABC transporter were found in response to both the imidacloprid and chlorpyrifos treatment. Trypsin had the third-largest number of DEUs after insecticide treatment for 36 h, 32 DEUs with 15 upregulated and 17 downregulated were observed in response to imidacloprid treatment, while 20 DEUs with eight upregulated and 12 downregulated were observed in the chlorpyrifos treatment. Additionally, the same 20 DEUs related to ABC transporter were found in both the imidacloprid and chlorpyrifos treatment. For the metabolism enzyme-related unigenes, such as cytochrome P450, GST, and carboxylesterase, there were two, two, and 11 cytochrome P450 related DEUs observed in response to chlorpyrifos treatment for 15 min, 3 h, and 36 h, respectively, and the same two DEUs with one upregulated and one downregulated were found at all treatment times. There were two, seven, and 11 cytochrome P450 related DEUs observed in response to imidacloprid treatment for 15 min, 3 h, and 36 h, respectively, with one identical upregulated DEUs found in all time durations, and one identical downregulated DEU found in the 15 min and 36 h groups, as well as the same three upregulated DEUs found after 3 h and 36 h of treatment. DEUs related to GST were only observed after 36 h of treatment, and the same downregulated DEU was found in response to both imidacloprid and chlorpyrifos treatment. DEUs related to carboxylesterase were only observed after 36 h of treatment, and the same six downregulated DEUs were observed in response to imidacloprid and chlorpyrifos treatment (Supplementary File S12). Finally, the expressed levels of one ABC transporter, one GST, one esterase, three cytochrome P450, two UGT, three trypsin related DEUs were chosen for qRT-PCR validation. Nine of sixty-six comparisons of RNA-Seq and qRT-PCR results do not agree, but all of them appear in the treatment of 15 min ( Figure 6).

Discussion
Insect resistance to pesticides is a complex impediment to agricultural production. Understanding how insects develop resistance to insecticides and the insecticide tolerance mechanism involved would help reduce and delay this process in insects. Using next-generation technologies to reveal tolerance and analyze insecticide-related genes by transcriptome profiles not only makes up for the gaps in previous studies but also provides us with new thoughts regarding insecticide resistance in S. avenae.
In this study, the number of DEUs increased as the insecticides' treatment time increased, and the KEGG and GO pathways were more abundant. Furthermore, some insecticide-related genes were only differentially expressed after treatment with insecticides for 36 h, including glutathione stransferase, carboxylesterase, acetylcholinesterase, acetylcholine receptor, chloride channel, and superoxide dismutase. Some genes may have been up-or downregulated when grain aphids were under insecticide pressure for 15 min and 3 h, but not significantly. After 36 h of insecticide treatment, they may have significantly up-or downregulated. In B. odoriphaga, when samples were treated with

Discussion
Insect resistance to pesticides is a complex impediment to agricultural production. Understanding how insects develop resistance to insecticides and the insecticide tolerance mechanism involved would help reduce and delay this process in insects. Using next-generation technologies to reveal tolerance and analyze insecticide-related genes by transcriptome profiles not only makes up for the gaps in previous studies but also provides us with new thoughts regarding insecticide resistance in S. avenae.
In this study, the number of DEUs increased as the insecticides' treatment time increased, and the KEGG and GO pathways were more abundant. Furthermore, some insecticide-related genes were only differentially expressed after treatment with insecticides for 36 h, including glutathione s-transferase, carboxylesterase, acetylcholinesterase, acetylcholine receptor, chloride channel, and superoxide dismutase. Some genes may have been up-or downregulated when grain aphids were under insecticide pressure for 15 min and 3 h, but not significantly. After 36 h of insecticide treatment, they may have significantly up-or downregulated. In B. odoriphaga, when samples were treated with chlorpyrifos and clothianidin for 6 h and 48 h, the number of DEUs related to insecticide tolerance was not much different, and the majority of insecticide tolerance-related unigenes responded in 6 h; however, there were still some insecticide tolerance-related unigenes that responded in 48 h [6], which was similar to the results of the current study. Therefore, the results of this study may indicate that exposure to insecticides for 3 h does not lead to a great increase in tolerance-related unigenes. After 36 h of insecticide treatment, many more insecticide tolerance-related unigenes were differentially regulated, including those that responded to short-term exposure. These results indicated that the treatment duration of insecticides has a greater impact on DEUs related to insecticide tolerance than the type of insecticide.
Cytochrome P450 is an enzyme that has a variety of metabolic functions, and increased cytochrome P450 mediated drug metabolism is an important detoxification mechanism for insects [22]. The overexpression of the P450 monooxygenase enzyme is the most common mechanism of imidacloprid and chlorpyrifos resistance [23][24][25]. In aphids, P450 monooxygenase enzyme also plays an important role in insecticide detoxification and resistance [26][27][28]. In the present study, samples treated with chlorpyrifos showed five upregulated P450 unigenes, as well as a fold change in imidacloprid-treated samples. Four of the five unigenes were annotated as the CYP4Csubfamily gene, and one was annotated as the CYP6A subfamily gene. Unigene DN67665_c1_g4 showed upregulation with time, and unigene DN68255_c7_g3 showed the highest upregulated fold change (3.29). Furthermore, the five aforementioned upregulated CYP unigenes belonged to the CYP 4 and 6 families, and those two family genes have been implicated in insecticide resistance more often than any other P450 family [22,29]. These genes are unique to insects and play important roles in the metabolism and detoxification of pesticides [22], and it has been suggested that overexpression of unigenes in these two families was involved in insecticide tolerance and detoxification in S. avenae. In B. odoriphaga, CYP6FV12 showed different fold changes in different life stages when exposed to imidacloprid and was confirmed to be related to B. odoriphaga resistance to imidacloprid [30], while CYP6CM1vQ was confirmed to be associated with a high level of imidacloprid resistance in Bemisia tabaci [31]. Four other upregulated unigenes were only differentially expressed following treatment with imidacloprid for 3 h, and all of four of these were annotated as the CYP6A subfamily gene, which indicated that this subfamily gene may be important to S. avenae detoxification to imidacloprid. Therefore, although our results demonstrated that overexpression of CYP 4 and 6 family genes associated with the detoxification of imidacloprid and chlorpyrifos in S. avenae, whether these P450s can metabolize imidacloprid and chlorpyrifos needs further research.
ABC transporters play an important role in the detoxification process phase III, which can transport the polar compounds or conjugates out of the cell [32]. ABC transporters have been associated with imidacloprid and chlorpyrifos resistance in insects [33,34]. Eight ABCB/C/D/G subfamily transporter genes in imidacloprid and chlorpyrifos resistant strains of Laodelphax striatellus were significant upregulated compared with a susceptible strain [35], these results suggest that ABC transporters might be involved in resistance to multiple insecticides in L. striatellus. Two out the five ABC transporter genes analyzed in Anopheles gambiae were downregulated after the 48h exposure of permethrin [36]. In our study, six of the nine up-regulated ABC transporters belonged to C and G subfamilies. ABC transporter possessed a larger group of DEUs induced by imidacloprid and chlorpyrifos than the groups of P450, GST, and carboxylesterase, and more than two of three ABC transporter unigenes are downregulated after pesticides treatments. Thus, ABC transporter may play an important role in the detoxification process and insecticide tolerance of S. avenae.
In addition to insecticide detoxification, target site sensitivity and decreased penetration are important to insect tolerance to pesticides. In this study, no DEUs related to target site sensitivity were found, but dozens of penetration related DEUs were identified. Increased insecticides cuticular penetration, including cuticle thickening and alteration of cuticle composition, have previously been described [37]. In Culex pipiens pallens, cuticle protein played an important role in deltamethrin resistance [38], and CPLCG5 encoded a cuticle protein that participated in pyrethroid resistance by inducing rigidity and increasing the thickness of the cuticle [39]. There were five cuticle protein genes differentially expressed in deltamethrin-resistant C. pipiens pallens when compared with susceptible strains, with cuticle protein CP14.6 precursors found to be overexpressed in the deltamethrin-resistant strain. This may support the hypothesis that mosquitoes can protect themselves from insecticides by regulating cuticles, which finally leads to cuticular resistance [40]. In the present study, after exposure to insecticides, the largest group of DEUs was cuticle protein-related unigenes, and most of them were down regulated. This is also happening in C. pipiens pallens. In the 30 differentially expressed proteins identified by deltamethrin-resistant strain compared with the deltamethrin-susceptible strain, five out of 30 proteins are cuticle-related protein, and four out of five are downregulated [40]. Mevinphos resistance strains comparing with susceptible strains in P. xylostella, 12 out of 16 differentially expressed cuticle protein transcripts are downregulated [41]. Therefore, our study indicates that cuticle proteins may play an important role in metabolism or tolerance to insecticides by S. avenae, but how cuticle proteins are involved in the process of cuticle alterations, its alterations of cuticle structure or composition, and how to slow down the penetration of insecticides requires further research.
Trypsin-related genes accounted for a large group of DEUs in B. oriphaga under insecticide stress [6] and were found to be highly expressed in C. pipiens pallens deltamethrin-resistant strains [42]. Following exposure to triazophos, imidacloprid, chlorpyrifos, and abamectin, trypsin expression was upregulated in Sodatella furcifera [43]. In the present study, trypsin-related genes accounted for the largest group of upregulated unigenes under insecticide stress for 36 h, indicating that trypsin may be related to the response to stress induced by insecticides in S. avenae. In spirotetranmat-and thiamethoxam-resistant strains of A. gossypii, UGT was significantly upregulated relative to the susceptible strains [44,45]. In the present study, most of the UGTs were downregulated, while only one UGT was upregulated after exposure to insecticides, and all of the DEUs of UGT were assigned to KEGG as drug metabolism genes. In a female Spodoptera littoralis, exposure to a pheromone or plant odorant led to differential downregulation of the transcription levels of two UGTs specifically expressed in antennae [46]. Because UGTs played an important role in a variety of physiological and biochemical processes in insects, including detoxification of substrates (such as plant allelochemicals and insecticides) [47,48], our results indicated that UGTs may play a role in the tolerance and detoxification of insecticides; however, further study is needed to confirm these findings.
Sequencing by treatment with an insecticide to observe up-or downregulation of certain enzyme or receptor genes is the first step in understanding whether they are involved in pest resistance. Based on the analysis of unigenes that showed significant differences in the expression in response to these pesticides, it is concluded that the types and numbers of DEUs increased with the increased treatment time, while the differential unigene expression in response to different agents at the same time did not vary greatly. These findings may indicate that the production of S. avenae tolerance of insecticides does not occur via regulation by a single gene, but rather a result of joint regulation by multiple genes.

Conclusions
In this study, the adult transcriptome of S. avenae was sequenced, after which the unigenes database was assembled and this is the first time the annotation to different databases in S. avenae has occurred. The unigenes involved in responding to two insecticides, chlorpyrifos and imidacloprid, after different exposure times, were then identified and analyzed. The transcriptome assembly results provide a substantial contribution to the existing sequence resources for S. avenae. The analysis of DEUs responding to insecticides could provide a substantial foundation for research regarding the tolerance and detoxification mechanisms of S. avenae. The upregulated expression of cytochrome P450 genes may be important to pesticide detoxification in S. avenae. However, further investigation of the DEUs related to insecticide tolerance and detoxification is needed to determine if they can be used as molecular targets to explore novel approaches to control S. avenae.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/10/12/951/s1, File S1: R code for detect differential expression of unigenes; File S2: The primers to confirm unigenes and qRT-PCR; File S3: The amplification efficiencies of primers for qRT-PCR; File S4: Dose-response curves of S. avenae to chlorpyrifos and imidacloprid; File S5: Statistics of clean reads after RNA-Seq; File S6: Size distribution of unigenes; File S7: KEGG annotation and pathways of the S. avenae transcriptome; File S8: Cluster of orthologous groups functional classification of the S. avenae transcriptome; File S9: DEUs up-or downregulated in response to different treatments; File S10: GO enrichment analysis of DEUs; File S11: KEGG enrichment analysis of DEUs; File S12: Insecticide tolerance-related unigenes.