Transcriptome Profiling Revealed Potentially Critical Roles for Digestion and Defense-Related Genes in Insects’ Use of Resistant Host Plants: A Case Study with Sitobion Avenae

Using host plant resistance (HPR) in management of insect pests is often environmentally friendly and suitable for sustainable development of agricultural industries. However, this strategy can be limited by rapid evolution of insect populations that overcome HPR, for which the underlying molecular factors and mechanisms are not well understood. To address this issue, we analyzed transcriptomes of two distinct biotypes of the grain aphid, Sitobion avenae (Fabricius), on wheat and barley. This analysis revealed a large number of differentially expressed genes (DEGs) between biotypes 1 and 3 on wheat and barley. The majority of them were common DEGs occurring on both wheat and barley. GO and KEGG enrichment analyses for these common DEGs demonstrated significant expression divergence between both biotypes in genes associated with digestion and defense. Top defense-related common DEGs with the most significant expression changes included three peroxidases, two UGTs (UDP-glycosyltransferase), two cuticle proteins, one glutathione S-transferases (GST), one superoxide dismutase, and one esterase, suggesting their potentially critical roles in the divergence of S. avenae biotypes. A relatively high number of specific DEGs on wheat were identified for peroxidases (9) and P450s (8), indicating that phenolic compounds and hydroxamic acids may play key roles in resistance of wheat against S. avenae. Enrichment of specific DEGs on barley for P450s and ABC transporters suggested their key roles in this aphid’s detoxification against secondary metabolites (e.g., alkaloids) in barley. Our results can provide insights into the molecular factors and functions that explain biotype adaptation in insects and their use of resistant plants. This study also has significant implications for developing new resistant cultivars, developing strategies that limit rapid development of insect biotypes, and extending resistant crop cultivars’ durability and sustainability in integrated management programs.


Introduction
In the constantly changing agricultural landscape, effective management of various insect pests often presents huge challenges. Many pest management tactics, such as cultural, biological and chemical tactics, have been developed to meet these challenges. Among different strategies, deployment of host plant resistance can provide many benefits, such as economic savings, promotions of natural and biological control, and a decrease in potentially hazardous insecticide applications [1,2]. Thus, this pest management strategy developed from naturally evolved plant defenses is often environmentally friendly and suitable for sustainable development of agricultural industries, and it can be very valuable in the integrated management of insect pests for various agricultural crops [1,3]. However, some factors, such as efficacy, lack of resistant varieties and rapid evolution of insect populations, can limit the utilization of host plant resistance in insect pest management [2,4]. One of the most serious challenges to the use of resistant crops against various insect pests is the rapid development of new biotypes [5]. For example, a new biotype of the soybean aphid has developed in North America within 5 years after invasion, despite the genetic bottleneck during the invasion process [2]. Multiple biotypes have been constantly found in many phytophagous insects in response to different levels of host plant resistance, such as the rice brown plant hopper (Nilaparvata lugens), Hessian fly (Mayetiola destructor), and soybean aphid (Aphis glycines) [6][7][8].
In the context of host plant resistance, insect biotypes are defined as populations that are able to survive, reproduce on, and cause injury to a cultivated plant resistant to other populations of the same species [2,4]. So, biotypes can be determined by their unique response profiles (i.e., differential survival or fitness) on different plants [2,4]. Particularly, aphids are prone to adapt and overcome host plant resistance and develop different biotypes on variable plants. For example, fifteen sympatric biotypes of the pea aphid have been found on different host species, and they present gradual levels of genetic divergence [9,10]. A study has found at least 17 aphid species with multiple biotypes evolved [11]. This is probably because differential levels of resistance in variable plants present a huge pressure of natural selection that favors host-associated ecological divergence of various aphid populations [12,13]. Increasing evidence has shown that the differential levels of plant resistance to aphids can be attributed to plant secondary chemistry. For example, different levels of hydroxamic acids in cereal crops (e.g., wheat) can play a key role in resistance of crops against a wide range of insects, including aphids [14,15]. Facing harmful effects of plant secondary compounds, aphids have developed a diversity of detoxification systems, such as cytochrome P450 monooxygenases (P450s), glutathione S-transferases (GSTs) and esterases (ESTs) [16][17][18][19]. There are at least two possible means by which detoxification can lead to differential plant resistance and aphid biotype adaptation: (1) detoxification enzymes can be produced in differential ratios and amounts on different resistant plants; (2) gene mutations can occur in the catalytic sites of detoxification enzymes enabling a much more efficient and effective process of detoxification [2]. Activities of the above-mentioned defense-related enzymes in aphids have been extensively studied, but the underlying molecular variation has received little attention [20][21][22][23]. Since host plant-associated differentiation may often involve genes targeted by selection, identification of elevated levels of molecular differentiation between biotypes (or host-associated populations) is an important step towards fully understanding adaptive divergence, as well as the underlying mechanisms of host plant resistance [24]. Research in this respect can not only expand our understanding of insect-plant interactions, but also help to improve the efficacy, durability, and sustainability for the use of host plant resistance in management of insect pests [2]. So far, such research is still scarce, and molecular factors and functions underlying the use of resistant plants by different aphid biotypes (differential adaptation to different resistant plants) are not well understood.
The English grain aphid, Sitobion avenae (Fabricius) (Hemiptera: Aphididae), is a good model to address this issue. This aphid can feed and survive on many species of the Poaceae, including cereal crops and pasture grasses [3,[25][26][27]. These host plants exhibit a wide range of concentrations for secondary metabolites (e.g., phenolic compounds, hydroxamic acids and alkaloids), which are chemical defenses involved in resistance against aphids [28,29]. Our previous studies did find some genetic variation among S. avenae populations from different cultivated hosts and geographic areas, which could serve as a possible genetic reservoir for adaptation to differential host plant resistance [30][31][32][33]. In addition, based on their unique performance profiles on resistant wheat and Insects 2020, 11, 90 3 of 22 barley cultivars, multiple S. avenae biotypes were identified [4]. Significant genetic differentiation was detected among these biotypes (e.g., biotype 1 vs. biotype 3, F ST = 0.157) based on microsatellite data [33]. Despite this, molecular factors and functions underlying the divergence of biotypes in S. avenae are still little understood. In this study, molecular variation in two distinct S. avenae biotypes (i.e., biotypes 1 and 3) feeding on resistant host plants is determined by using high-throughput sequencing techniques. The objectives are to: (1) examine potential molecular factors underlying the adaptive divergence of S. avenae biotypes; (2) identify functions and genes involved in the use of resistant plants (i.e., adaptation to specific plants) by different S. avenae biotypes.

Aphid Sample Collection and Colony Establishment
In our previous study on the identification of S. avenae biotypes [4], we found that biotype 3 was characterized by an ability to overcome the resistance of barley (Hordeum vulgare L.) cultivars (e.g., Xiyin No.2), but not wheat (Triticum aestivum L.) cultivars (e.g., Aikang 58). On the contrary, biotype 1 was able to overcome the resistance of wheat cultivars (e.g., Aikang 58), but not barley cultivars (e.g., Xiyin No.2). In 2016, clones of biotypes 1 and 3 were sampled on wheat and barley, respectively [4]. The two biotypes of S. avenae were reared on the plant of origin (i.e., wheat or barley) under a temperature of 22 ± 1 • C, a relative humidity of 65 ± 5%, and a photoperiod of 16:8 (L:D) h. Prior to the following experiments, all aphid clones were maintained under common laboratory conditions for at least three generations for the purpose of minimizing confounding environmental effects.

Fitness Bioassays of Both S. Avenae Biotypes on Wheat and Barley
In order to determine their fitness on wheat and barley, new-born first instar nymphs of both S. avenae biotypes were transferred onto single two-leaf stage seedlings (one nymph per seedling) of Aikang 58 (i.e., wheat) and Xiyin No.2 (i.e., barley) planted in 200 mL plastic pots [6 cm in diameter, containing turfy soil mixed with vermiculite and perlite (4:3:1, v/v/v)]. Each plastic pot was well enclosed with a transparent plastic cylinder (6 cm in diameter, 15 cm in height) which had a terylene mesh top for ventilation. Test aphid individuals on both plants, they were kept in growth chambers with the above-mentioned conditions and monitored for molting and reproductive events daily. New-born aphid individuals were counted and then removed from the plant each day. This experiment was terminated on day 10 after the initiation of reproduction for each test aphid individual. Thirty replicates were conducted for each S. avenae biotype on each plant (i.e., wheat or barley). The 10 d fecundities (offspring accumulated in 10 days since the initiation of reproduction) for each biotype on both plants were tabulated and analyzed by using one-way ANOVA in SAS [34]. Post-hoc comparisons between treatments were conducted by using Tukey tests at α = 0.05 following significant ANOVA.

RNA Sampling and Sequencing
New-born individuals of biotypes 1 and 3 were transferred onto Aikang 58 and Xiyin No.2 seedlings (one nymph per seedling) at the two-leaf stage. They were kept under the aforementioned laboratory conditions until reaching the adult stage. After that, they were allowed to feed for additional 24 h, and sampled for RNA sequencing. Ten wingless aphid individuals were collected and put into a 1.5 mL RNase-free tube each time. All RNA-seq aphid samples were frozen immediately in liquid nitrogen and stored in a freezer at -80 • C until use in RNA extraction. There were three biological replications for each S. avenae biotype on each plant. The total RNA of each sample was extracted with the MiniBEST Universal RNA Extraction Kit (Takara Bio Inc., Dalian, China), and RNase-free DNase I (Takara Bio Inc., Dalian, China) was used to remove the potential genomic DNA contamination of samples. Following the manufacturers' instructions, the quality and quantity of RNA samples were estimated with a Bioanalyzer 2100 instrument (Agilent Technologies, CA, US) and NanoPhotometer ® spectrophotometer (IMPLEN, CA, US). The cDNA libraries for RNA samples were generated with the NEBNext ® UltraTM RNA Library Prep Kit for Illumina (NEB, Beverly, MA, US). The Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, US) was used to check the quality of these cDNA libraries. The paired-end DNA sequencing of all samples was conducted on an Illumina HiSeq 2500 platform (Illumina Inc., San Diego, CA) of Ovidson Gene Technology Co., Ltd., Beijing, China. Raw datasets of the above-mentioned RNA-seq samples have been submitted to the Sequence Read Archive (SRA) database of NCBI (National Center for Biotechnology Information) (NCBI accession ID: PRJNA575173).

Transcriptome Assembly and Annotation
Through filtering out low quality reads (more than 50% of nucleotides with Qphred ≤ 20), and reads with adapters and ambiguous "N" nucleotides >10% in raw datasets, clean datasets were obtained. At the same time, Q20, Q30 and GC-content of these clean datasets were calculated. We pooled clean reads of all samples and conducted the transcriptome assembly in TRINITY (v2.1.1) with default parameters (for details, see [35]). In order to minimize the redundancy effects, the software CD-HIT (v4.6.7) was used to process transcripts with 95% similarity [36]. The longest transcript of each gene was defined as a 'unigene' for functional annotation. Seven databases were used to annotate all unigenes: NR (NCBI non-redundant proteins, NCBI BLAST (basic local alignment search tool) 2.2.28+, e-value ≤ 1.0 × 10 −5 ), NT (NCBI non-redundant nucleotides, NCBI BLAST 2.2.28+, e-value ≤ 1.0 × 10 −5 ), Pfam (protein families, HMMER 3.0 package, hmmscan, e-value ≤ 0.01), KOG (eukaryotic orthologous groups, NCBI Blast 2.2.28+, e-value ≤ 0.001) and Swiss-Prot (a manually annotated and reviewed protein sequence database, NCBI BLAST 2.2.28+, e-value ≤ 1.0 × 10 −5 ). Further orthology analyses for unigenes were performed to determine their molecular functions by using the KEGG (Kyoto Encyclopedia of Gene and Genomes) automatic annotation server (KAAS, KEGG Automatic Annotation Server, e-value ≤ 1.0 × 10 −10 ) [1], and the resulting KO (KEGG Orthologs) terms were compared between the two biotypes. Gene ontology (GO) terms for each unigene were further processed in Blast2GO v2.5 with a threshold e-value of ≤ 1.0 × 10 −5 [37].

Analysis of Gene Expression Patterns
The expression levels of each unigene were determined by RSEM v 1.2.3 (bowtie2 default parameters, mismatch = 0) through mapping each clean read to the assembled transcriptome [38]. Based on the obtained matrix of raw counts for each unigene in each sample, the differential expression of genes (DEGs) between two S. avenae biotypes on each plant (i.e., wheat or barley) was analyzed with DESeq2 (v 1.10.1) by using default parameters [39]. We identified 3392 and 2246 DEGs between biotypes 1 and 3 on wheat and barley, respectively. DEGs co-existed in aphids feeding on wheat and barley were considered to be common DEGs (i.e., 1528)-otherwise, they were considered wheat-specific DEGs (i.e., 1864) or barley-specific DEGs (i.e., 718). The significance of enrichment of differential expression genes in KEGG pathways was tested by using the KOBAS software [40]. GO term enrichment analysis for these DEGs were analyzed and visualized by using the online tools of GO analysis (https://www.omicshare.com/tools/Home/Soft/gogsea).

Quantitative Real-Time PCR (qRT-PCR)
In order to verify the results of the above-mentioned RNA-seq analysis, the expression of 12 selected genes was examined by using the qRT-PCR (quantitative reverse transcription-polymerase chain reaction) method. The gene NADH (nicotinamide adenine dinucleotide dehydrogenase, c93747_g4) was chosen as the reference gene because of its consistent expression in our study [M (stability measure value) = 0.177, SV (stability value) = 0.039]. Gene-specific primers were designed by using Beacon Designer version 8.0 (Premier Biosoft, Palo Alto, CA). All the primers used are listed in Table S1. PCR reactions were prepared with the SYBR Premix Ex Taq™ Kit (Takara, Dalian, China) according to the manufacturer's instructions. All qRT-PCR reactions were performed on the Roche LightCycler ® 480 II system (Roche Diagnostics Ltd., Rotkreuz, Switzerland), with three biological replicates and two technical replicates for each gene. The qPCR cycling parameters were as follows: 95 • C for 5 min, followed by 40 cycles of 95 • C for 10 s and 60 • C for 30 s. To confirm the homogeneity of the PCR products, reaction stages for establishment of melt curves (95 • C for 5 s, 65 • C for 1min, and 95 • C for 5 s) were also included. The relative expression level of each gene was calculated by using the 2 −∆∆Ct method [41]. The relative expression levels and FPKM values of each gene were compared between the two biotypes by using the student t test in the software SPSS Statistics 23 (IBM SPSS Statistics Inc., Chicago, IL, USA).
Insects 2020, 11, x FOR PEER REVIEW 5 of 23 LightCycler ® 480 II system (Roche Diagnostics Ltd., Rotkreuz, Switzerland), with three biological replicates and two technical replicates for each gene. The qPCR cycling parameters were as follows: 95 °C for 5 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s. To confirm the homogeneity of the PCR products, reaction stages for establishment of melt curves (95 °C for 5 s, 65 °C for 1min, and 95 °C for 5 s) were also included. The relative expression level of each gene was calculated by using the 2 −ΔΔCt method [41]. The relative expression levels and FPKM values of each gene were compared between the two biotypes by using the student t test in the software SPSS Statistics 23 (IBM SPSS Statistics Inc., Chicago, IL, USA).

De Novo Assembly and Transcriptome Annotation
In total, 343,816,452 raw reads and 337,188,919 clean reads were obtained in this study ( Table 1). The Q20 of all samples was >98.66%, and the Q30 was >96.18%, showing the high quality of the RNA-Seq analysis. Each sample library was mapped back to the full assembly, with overall alignment rates of 71.67%-73.23%. The de novo assembly yielded 143,058 unigenes. The median length of these unigenes was 358 bp, with a N50 (length N for which 50% of all bases in the assembly are located in a unigene of length < N) of 1012 bp. As shown in Figure 2A, the unigenes of length < 300 bp and >500 bp accounted for 39.70% and 34.70% of the total, respectively, whereas 25.60% of all the unigenes fell between 300 and 500 bp.

De Novo Assembly and Transcriptome Annotation
In total, 343,816,452 raw reads and 337,188,919 clean reads were obtained in this study ( Table 1). The Q20 of all samples was >98.66%, and the Q30 was >96.18%, showing the high quality of the RNA-Seq analysis. Each sample library was mapped back to the full assembly, with overall alignment rates of 71.67%-73.23%. The de novo assembly yielded 143,058 unigenes. The median length of these unigenes was 358 bp, with a N50 (length N for which 50% of all bases in the assembly are located in a unigene of length < N) of 1012 bp. As shown in Figure 2A, the unigenes of length < 300 bp and >500 bp accounted for 39.70% and 34.70% of the total, respectively, whereas 25.60% of all the unigenes fell between 300 and 500 bp. Note: Q20 and Q30, percentage of bases with quality over a Phred score of 20 (i.e., an error rate of 1%) and 30 (i.e., an error rate of 0.1%), respectively; GC, proportion of G and C bases; N50, length N for which 50% of all bases in the assembly are located in a unigene of length < N.

Identification of Common and Specific DEGs (Differentially Expressed Genes) and Enrichment Analyses
On wheat, 3392 DEGs (adjusted p values <0.05, fold change values >2) were detected between S. avenae biotypes 1 and 3 ( Figure 3). For these DEGs, 1513 of them were up-regulated, and 1879 were down-regulated in biotype 3 compared to biotype 1 (biotype 3 vs. biotype 1). On barley, 1048 up-and 1198 down-regulated genes (2246 DEGs in total) were detected between the two biotypes (biotype 3 vs. biotype 1). In addition, 1528 DEGs (referred to as common DEGs hereafter) occurred in the two biotypes feeding on both wheat and barley ( Figure 4A), and they included 695 up-regulated genes and 833 down-regulated genes (biotype 3 vs. biotype 1, Figure 4B). One thousand eight hundred sixty-four and 718 DEGs (referred to as specific DEGs hereafter) occurred only in the two biotypes feeding on wheat and barley, respectively ( Figure 4A). Of the 1864 specific DEGs on wheat, 818 and 1046 genes were up-regulated and down-regulated, respectively (biotype 3 vs. biotype 1) ( Figure 4B). Of the 718 specific DEGs on barley, 353 and 365 genes were up-regulated and down-regulated, respectively (biotype 3 vs. biotype 1) ( Figure 4B).
For specific DEGs on barley, 11 genes related to digestion were identified, including one serine protease, three lipases, one phospholipase, five maltases, and one serine carboxypeptidase (Figure 7). Among them, all three lipases, one maltase, the serine protease and the serine carboxypeptidase were up-regulated in biotype 3 compared with biotype 1. Thirteen defense-related specific DEGs on barley were also identified, including four P450s, four UGTs, one ABC transporter, one POD, two CPs, and one trehalose transporter. Of them, all four P450s, two UGTs, the ABC transporter, and one trehalose transporter were up-regulated in biotype 3 compared with biotype 1. Among these, the top defenserelated specific DEGs (biotype 3 vs. biotype 1) on barley with the most significant changes included two UGTs [i.e., c80170_g1 and c82523_g1 with log2(fold change) values being 5.39 and −1.53, respectively], two P450s [i.e., c70940_g1 and c78468_g1 with log2(fold change) values being 2.71 and 1.58, respectively], and one CP [i.e., c92386_g6, log2(fold change): 2.12]. Figure 6. A heatmap of specific DEGs related to digestion and defense in two Sitobion avenae biotypes on wheat (genes with expression higher and lower than the mean are indicated by blue and yellow, respectively).
For specific DEGs on barley, 11 genes related to digestion were identified, including one serine protease, three lipases, one phospholipase, five maltases, and one serine carboxypeptidase (Figure 7). Among them, all three lipases, one maltase, the serine protease and the serine carboxypeptidase were up-regulated in biotype 3 compared with biotype 1. Thirteen defense-related specific DEGs on barley were also identified, including four P450s, four UGTs, one ABC transporter, one POD, two CPs, and one trehalose transporter. Of them, all four P450s, two UGTs, the ABC transporter, and one trehalose transporter were up-regulated in biotype 3 compared with biotype 1. Among these, the top defense-related specific DEGs (biotype 3 vs. biotype 1) on barley with the most significant changes included two UGTs [i.e., c80170_g1 and c82523_g1 with log2(fold change) values being 5.39 and −1.

Divergence of Aphid Biotypes
Host plant resistance has been used in many integrated management programs for insect pests. However, its use is often limited by the rapid development of various insect biotypes. Molecular researches involving insect biotypes may provide insights into the underlying molecular factors and mechanisms for evolution of biotypes [13,42]. In this study, we used the English grain aphid (Sitobion avenae) as a model to address these issues. This aphid was found to have evolved multiple biotypes [4,43]. Among them, biotypes 1 and 3 were showed to have developed a certain degree of specialization on wheat and barley, respectively [4]. The two S. avenae biotypes showed differential performance on wheat (Aikang 58) and barley (Xinyin No.2) in this study, confirming their adaptive divergence on wheat and barley reported in our previous study [4]. This is also consistent with the finding that some S. avenae clones showed different degrees of specialization on barley and wheat in Figure 8. Validation of RNA-Seq results with qRT-PCR of 12 unigenes (A), and Pearson's correlation between fold changes from both qRT-PCR and RNA-Seq (B) (c84549_g2, cytochrome P450 4C1; c94527_g1, UDP-glucuronosyltransferase 2B7; c92512_g3, esterase E4; c74919_g1, zinc transporter ZIP1-like; c28877_g1, heat shock protein 83; c75571_g2, cysteine protease ATG4B; c70940_g1, cytochrome P450 6a13; c86074_g2, cytochrome P450 307a1-like; c80170_g1, UDP-glucuronosyltransferase 2C1; c89581_g2, ABC transporter G family member 20; c75250_g2, trehalose transporter Tret1-like; c79176_g2, serine protease 44-like; the first and last six genes were randomly selected from specific DEGs on wheat and barley, respectively; * p < 0.05; ** p < 0.01; *** p < 0.001).

Divergence of Aphid Biotypes
Host plant resistance has been used in many integrated management programs for insect pests. However, its use is often limited by the rapid development of various insect biotypes. Molecular researches involving insect biotypes may provide insights into the underlying molecular factors and mechanisms for evolution of biotypes [13,42]. In this study, we used the English grain aphid (Sitobion avenae) as a model to address these issues. This aphid was found to have evolved multiple biotypes [4,43]. Among them, biotypes 1 and 3 were showed to have developed a certain degree of specialization on wheat and barley, respectively [4]. The two S. avenae biotypes showed differential performance on wheat (Aikang 58) and barley (Xinyin No.2) in this study, confirming their adaptive divergence on wheat and barley reported in our previous study [4]. This is also consistent with the finding that some S. avenae clones showed different degrees of specialization on barley and wheat in previous studies [3,30]. We further tested the two S. avenae biotypes on both wheat and barley by using transcriptome profiling techniques. A large number of DEGs were identified between the two biotypes when feeding on wheat (i.e., 3392) or barley (i.e., 2246). Interestingly, the majority of them (i.e., 1528) were common DEGs (occurring on both wheat and barley). The enriched GO terms for these common DEGs demonstrated that there had been expression divergence between both biotypes in genes associated with digestion (e.g., peptidase, serine-type peptidase, and serine hydrolase) and defense (e.g., UDP-glycosyltransferase). Similarly, in KEGG enrichment analyses, some digestion-related (e.g., bile secretion) and detoxification-related pathways (e.g., drug and xenobiotics metabolism) were also significantly enriched. Based on these analyses, we identified 16 common DEGs related to digestion, and 40 to defense that could be involved in the biotype divergence process in S. avenae.
The digestion-related common DEGs encoded five lipases, three phospholipases, two carbohydrases, one maltase, one galactosidase, one mannosidase, and three serine proteases. It Tmakes sense since many of them (e.g., maltase, galactosidase, mannosidase, lipase and phospholipase) are essential digestive enzymes for glycolysis or esterlysis, and critical to the development and reproduction of insects [44][45][46][47][48][49]. Differential expression of a significant number of digestion-related genes was also identified among different biotypes in the whitefly Bemisia tabaci and fruit fly Drosophila patchea [50,51]. Thus, the differential expression of these digestive enzymes may have significant implications in the adaptation of insect biotypes on different host plants with significant nutritional variation [52][53][54][55]. In addition to their nutritional implications, some proteases, such as the above-mentioned serine proteases, may also be involved in defensive adaptation against proteinase inhibitors of host plants for insects [56,57]. Indeed, proteinase inhibitors of wheat have been proved to have anti-metabolic effects on S. avenae midgut proteases [58].
In addition to serine proteases, 40 other defense-related common DEGs were also identified between the two S. avenae biotypes. Among these, the top defense-related common DEGs (biotype 3 vs. biotype 1) with the most significant changes included three peroxidases (PODs), two UDP-glycosyltransferases (UGTs), two cuticle proteins (CPs), one glutathione S-transferase (GST), one superoxide dismutase (SOD), and one esterase (EST) ( Table 5). Previous studies have shown that defensive enzymes can be involved in the differentiation of pea aphid biotypes [13,55]. In addition, differences in the expression of genes related to detoxification were also detected between Subpsaltria yangi populations specialized on different host plants with variable compositions and concentrations of secondary plant compounds [59,60]. Therefore, all the above-mentioned defense-related common DEGs may play important roles in the divergence of S. avenae biotypes. Based on the 12 selected DEGs (Figure 8), nucleotide differences (p-distance = 0.006) showed little genetic divergence between the two biotypes. However, significant genetic differentiation between the two biotypes was detected based on microsatellite data in our previous study [33]. This makes sense since microsatellites are usually located in neutral regions, and the 12 selected DEGs can occur in conservative regions. In addition, in this study, the 12 selected DEGs showed large differences in expression between the two biotypes, indicating the importance of gene expression and regulatory components in the divergence of S. avenae biotypes. Further studies on specific functional roles of these defense-related DEGs and their regulatory components in S. avenae will increase our understanding of the phenomenon that biotype 1 had much higher proportions and much wider geographic distributions than biotypes 3 [4], as well as the process of biotype evolution in this aphid.

Molecular Factors Underlying Aphids' Use of Resistant Host Plants
Understanding the molecular differences among these apparently specialized aphid populations can have significant implications for deployment of host plant resistance in management of various pest aphids. In this study, S. avenae biotype 1 showed higher fecundity on wheat than on barley, and the opposite was true for biotype 3, suggesting that biotypes 1 and 3 could overcome the resistance of wheat and barley, respectively. We compared the transcriptomes of both S. avenae biotypes on the two plants and found that 1864 and 718 DEGs between both biotypes occurred only on wheat and barley, respectively (referred to as specific DEGs). The enriched GO terms for specific DEGs demonstrated that there were significant expression differences in genes associated with defense on wheat (e.g., those for activities of peroxidases, detoxification, and responses to toxic substance) and on barley (e.g., activities of oxidoreductases and gamma-glutamyltransferases). Similarly, in KEGG enrichment analyses of specific DEGs on barley, detoxification-related pathways (e.g., drug metabolism-cytochrome P450, metabolism of xenobiotics by cytochrome P450, drug metabolism-other enzymes, and ABC transporters) were also significantly enriched.
Our further screening analyses for defensive genes revealed 36 specific DEGs on wheat. Among them, a relatively high number of DEGs were identified for PODs (9), P450s (8), UGTs (4), and CPs (4). On barley, only 13 defensive DEGs were found to be specific, and they included four P450s, four UGTs, two CPs, one ABC transporter, one POD, and etc. Thus, wheat and barley with different levels of resistance to either test biotype induced specific expression of many defensive genes in S. avenae. For example, wheat induced expression changes for nine specific PODs, whereas only one POD was specifically induced by barley. The specific expression of defense-related genes in the two S. avenae biotypes may reflect differential compositions and contents of secondary metabolites (e.g., hydroxamic acids, alkaloids, phenolic compounds, and proteinase inhibitors) between wheat and barley reported in past studies [61][62][63][64][65][66][67]. Indeed, POD activities in S. avenae were found to be positively correlated with concentrations of diet phenolic compounds (e.g., flavonoids and tanins) [16,23]. In this study, three DEGs of PODs did show the most significant upregulation [log2 (fold change) for c81908_g5: 11.13; for c81908_g4: 9.34; for c81908_g1: 9.32] in S. avenae biotype 1 on wheat compared to biotype 3. The activity of P450s in S. avenae was associated with the concentration of plant hydroxamic acids [21,68,69]. High number of specific DEGs on wheat for PODs and P450s, as well as their high expression changes, may indicate that phenolic compounds and hydroxamic acids play critical roles in resistance of wheat against S. avenae. For barley, high concentrations of indole alkaloids (e.g., gramine) have often been found to be associated with host plant resistance to aphids [70][71][72][73][74]. In this study, GO and KEGG enrichment analyses showed that P450s and ABC transporters could be critical in defense of S. avenae on barley. Thus, P450s and ABC transporters can be key specific DEGs in metabolism of alkaloids in S. avenae on barley.
In addition to conventional detoxification DEGs, some digestive DEGs may also be important in S. avenae's use of resistant plants. For example, as one of the critical digestive enzymes responsible for cellulose degradation and glucose production, β-glucosidase was also found to be involved in aphids' metabolism of the plant secondary metabolite DIMBOA, which occurred frequently in wheat, instead of barley [75][76][77][78][79]. This explains the finding that the expression changes of β-glucosidase were specifically induced in S. avenae feeding on wheat in our study. Some proteases, like serine proteases and cysteine proteases, have both digestive and defensive functions in insects [56,57]. In this study, some cysteine proteases showed specific expression in S. avene on wheat, whereas some serine proteases were found to be specifically induced by barley. This suggests that serine and cysteine protease inhibitors can play important roles in resistance to aphids on barley and wheat, respectively. Further studies are needed to determine and verify the functional roles of all the putative key specific DEGs (digestive and/or defensive) in S. avenae against corresponding secondary compounds in wheat and barley.

Conclusions
In summary, this RNA-seq analysis provided a solid foundation toward improving our understanding of the molecular factors and functions underlying S. avenae's adaptation to different resistant host plants and its use of resistant hosts. Our study revealed a large number of digestive and defensive DEGs between S. avenae biotypes on both wheat and barley, suggesting that digestive and defensive enzymes could facilitate adaptations of S. avenae populations to changes in host plant chemistry. Different numbers of defensive DEGs were specifically induced by wheat and barley, reflecting differential compositions and contents of secondary metabolites (e.g., hydroxamic acids, alkaloids, phenolic compounds, and proteinase inhibitors) between both plants. A number of key defensive DEGs in S. avenae were found to be potentially involved in the metabolism of the above-mentioned plant secondary metabolites, and they can play a critical role in S. avenae's use of resistant host plants. Detailed studies are needed to determine the specific differences in the secondary metabolism of test plants, and exact functional roles of the putative key defensive DEGs in aphids' successful use of resistant plants. In addition, some salivary and odorant-binding proteins (c62817_g1 and c71932_g1) were also detected in common DEGs between the two S. avenae biotypes, suggesting the significance of these proteins in the differentiation of biotypes in aphids [13,24,42,55]. Further research with these putative key genes will help us to understand the molecular factors and mechanisms that explain virulence and biotype adaptation in insects. It is also a critical step in developing strategies that limit rapid development of insect biotypes and extending resistant crop cultivars' durability and sustainability in integrated management programs.