Transcriptome response comparison between vector and non-vector aphids after feeding on virus-infected wheat plants

Background Plant viruses maintain intricate interactions with their vector and non-vector insects and can impact the fitness of insects. However, the details of their molecular and cellular mechanisms have not been studied well. We compared the transcriptome-level responses in vector and non-vector aphids (Schizaphis graminum and Rhopalosiphum padi, respectively) after feeding on wheat plants with viral infections (Barley Yellow Dwarf Virus (BYDV) and Wheat dwarf virus (WDV), respectively). We conducted differentially expressed gene (DEG) annotation analyses and observed DEGs related to immune pathway, growth, development, and reproduction. And we conducted cloning and bioinformatic analyses of the key DEG involved in immune. Results For all differentially expressed gene analyses, the numbers of DEGs related to immune, growth, development, reproduction and cuticle were higher in vector aphids than in non-vector aphids. STAT5B (signal transducer and activator of transcription 5B), which is involved in the JAK-STAT pathway, was upregulated in R. padi exposed to WDV. The cloning and bioinformatic results indicated that the RpSTAT5B sequence contains a 2082 bp ORF encoding 693 amino acids. The protein molecular weight is 79.1 kD and pI is 8.13. Analysis indicated that RpSTAT5B is a non-transmembrane protein and a non-secreted protein. Homology and evolutionary analysis indicated that RpSTAT5B was closely related to R. maidis. Conclusions Unigene expression analysis showed that the total number of differentially expressed genes (DEGs) in the vector aphids was higher than that in the non-vector aphids. Functional enrichment analysis showed that the DEGs related to immunity, growth and reproduction in vector aphids were higher than those in non-vector aphids, and the differentially expressed genes related to immune were up-regulated. This study provides a basis for the evaluation of the response mechanisms of vector/non-vector insects to plant viruses.

secreted protein. Homology and evolutionary analysis indicated that RpSTAT5B was closely related to Rhopalosiphum maidis. Conclusions Unigene expression analysis showed that the total number of differentially expressed genes (DEGs) in the vector aphids was larger than that in the non-vector aphids.
Functional enrichment analysis showed that the DEGs related to immunity, growth and reproduction in vector aphids were larger than those in non-vector aphids, and the differentially expressed genes related to immune were up-regulated. This study provides a basis for the evaluation of the response mechanisms of vector/non-vector insects to plant viruses.

Background
Schizaphis graminum and Rhopalosiphum padi are serious pests of cereal crops, particularly wheat that may cause harm to plants by feeding on them and by transmitting the Barley yellow dwarf virus (BYDV) (Luteoviridae: Luteovirus) in a persistent and non-proliferative manner [1][2][3][4]. That causes one of the most economically important viral diseases of cereal plants [5,6]. In China, there were 4 isolated strains be identi ed as GPV, GAV, PAV and RMV and the BYDV-GAV is the major isolate [7]. Wheat dwarf virus (WDV) (Geminiviridae: Mastrevirus) is another serious virus of wheat in China that causes signi cant losses and is mainly transmitted by the leafhopper Psammotettix alienus (Dahlbom) in a persistent and nonproliferative manner [8,9]. S. graminum is a vector of BYDV-GAV and is a non-vector of WDV, while R. padi is a non-vector of both BYDV-GAV and WDV.
The aphids (S. graminum and R. padi) and viruses (BYDV-GAV and WDV) used in our study are important insect pests and viral pathogens that often occur together on the same host plant in agro-ecosystems. The plant-arthropod-pathogen interactions are complex, and some studies have shown that plant viruses can produce favorable or unfavorable effects on vector/or non-vector insects [10,11]. Exploring the effect of viruses on aphids can provide a better understanding of these three-way interactions and improve integrated pest management.
Previous studies showed that plant virus infection of its host plant could result in either bene cial or adverse effects on its vector via direct feeding or other indirect ways [12,13]. Transcriptional induction or repression is an important mechanism for insects to regulate their innate immune responses [14]. There are a few previous reports that have studied the global transcription pro les of insect vectors fed on virusinfected plants. Brault [15] conducted the rst large-scale analysis of aphid gene regulation following virus acquisition. Southern rice black-streaked dwarf virus (SRBSDV) infection activated the immune regulatory systems of the white-backed planthopper (WBPH) [16]. Transcriptome analysis of the thrips response to Tomato spotted wilt virus (TSWV) infection showed that the diversity of the innate immunerelated transcripts in response to viral infection was most pronounced in the P1 stage [17]. Reproduction, embryo development and growth were associated with upregulated contigs in virus-exposed Frankliniella fusca adults [18].
Compared with vectors, there are fewer studies on non-vectors. Although less studied, there is some evidence that plant viruses have an impact on non-vectors. The nymph survival rate was decreased and the longevity of female adults was shortened, while egg hatchability increased in response to non-vector brown planthoppers feeding on SRBSDV-infected plants [19]. In contrast, rice plants infected by rice black-streaked dwarf virus (RBSDV) improved the ecological tness of the brown planthopper [20]. In addition, TSWV infection reduced the fecundity and longevity of B. tabaci and increased the developmental time of Bemisia tabaci [21]. Additionally, feeding on plants inoculated with TSWV enhanced the developmental and oviposition rates of spider mites [22].
Pathogenic infections may affect the physiology of the host insect and even cause death [23]. In the process of long-term evolution, insects have formed defense mechanisms to achieve a dynamic balance between viruses and insects. Insects rely on their immune system to resist the invasion of pathogens, which is not conducive to the replication and transmission of the virus. One of the long-term goals of studying virus-insect interactions is to elucidate the molecular mechanism by which viruses affect the innate immune system of insects to avoid damage from the immune defense response. Therefore, understanding the insect's immune system is key to solving this problem. Immune-related signaling pathways that have been studied in insects include the RNAi, Toll, JAK/STAT, phagocytosis, apoptosis, proteolysis and JNK pathways [18,24,25].
Because the persistent plant viruses cannot break through the midgut barrier or salivary gland barrier of non-vector insects, they cannot be transmitted by the non-vector insects [26]. We previously have found that plant virus (e.g., BYDV-GAV and WDV) can be detected in non-vector aphids (e.g., R. padi) (unpublished data). Previous studies aslo have shown that plant viruses can improve the tness of nonvector insects (as well as vectors) [20][21][22][23]. We speculate that plant viruses in non-vector aphids can also affect their biochemical response. At present, the molecular and cytological mechanism of plant viruses to the vector and non-vector insects is not clear, so BYDV-GAV and WDV were used as the test viruses, R. padi and S. graminum were used as the test insects to compare the transcriptome changes of the vector and non-vector insects to the plant viruses, so as to nd out the differences between the plant viruses affecting the gene expression of the vector and non-vector insects. This study will lay the foundation for the molecular mechanism of plant viruses to improve the tness of non-mediator insects.

Sequence assembly and annotation
To explore the transcriptome response to viral infection by vector and non-vector aphids, we performed RNA-Seq analysis on adult S. graminum and R. padi exposed to BYDV-infected, WDV-infected and uninfected wheat plants. The S. graminum and R. padi cDNA libraries were sequenced, which generated 72,092 and 68,996 unigenes, respectively ( Table 1). The N50 sizes were 1,481 kb and 1,497 kb, respectively. The Q30 of each transcriptome library was above 85.74% and 85.62% for S. graminum and R. padi, respectively.
Nr (non-redundant protein sequence database) homologous species distribution revealed that 71.64% of S. graminum sequences matched with Acyrthosiphon pisum ( Fig. 1 A), and 67.81% of R. padi sequences matched with A. pisum ( Fig. 1 B). For the Gene ontology assignment, 7263 and 8164 unigenes in S. graminum and R. padi could be annotated. For S. graminum, 6200 unigenes were associated with biological process, 3080 with cellular component, and 6344 with molecular function (Fig. S1 A). In R. padi, 7104 unigenes were associated with biological process, 3597 with cellular component, and 7173 with molecular function (Fig. S1 B). These sequences have been submitted to the NCBI Sequence Read Archive (SRA), and the accession number is PRJNA490258. DEGs in S. graminum and R. padi, as the vector/non-vector, in response to feeding on virus-infected wheat plants For S. graminum fed on BYDV-infected wheat, 1525 DEGs were identi ed, with 693 upregulated and 832 downregulated. There were fewer DEGs in the S. graminum fed on WDV-infected wheat (494), and the majority were downregulated (73.28%). There were 589 DEGs in R. padi exposed to BYDV-infected wheat, with 247 upregulated and 342 downregulated. A total of 343 DEGs were identi ed in R. padi exposed to WDV-infected wheat, with 75.8% downregulated. The change ratios of most DEGs were between 2 -5 and 2 5 (Fig. 2). Gene ontology (GO) and KEGG enrichment analyses were conducted for all DEGs.
Distribution and enrichment of DEGs in vector and nonvector aphids into gene ontologies (GOs) Gene ontology (GO) analysis of all the DEGs in S. graminum and R. padi was conducted to classify the sequences into the biological process, molecular function and cellular component categories ( Fig. S2-Fig. S5). The most enriched category in both vector and non-vector aphids was metabolic process (Fig.  S2-Fig. S5). Among all comparisons, the differentially expressed genes contained a high proportion of genes involved in development, metabolism, reproduction and growth, and the number in vector aphids was higher (29 DEGs in Sg-BYDV) than that in non-vector aphids (13,13, and 3 DEGs in Sg-WDV, Rp-BYDV, and Rp-WDV, respectively) ( Fig. 3 and Table S5). In addition, genes associated with the cuticle (structural constituent of the cuticle and structural constituent of the chitin-based cuticle) were also differentially expressed, and there were 16 DEGs in the vector treatment (Sg-BYDV), while there were fewer in the non-vector treatments (with 2 in Sg-WDV, 5 in Rp-BYDV, and 2 in Rp-WDV). In Sg-BYDV, there were 15 DEGs associated with the cuticle, and the majority of them were downregulated (14 downregulated) ( Fig. 3 and Table 2). Several genes involved in the cytoskeleton were also differentially expressed, with 6, 1, and 2 related DEGs in Sg-BYDV, Sg-WDV, and Rp-BYDV, respectively ( Fig. 3 and Table   3). The results of GO enrichment analysis showed that the number of categories of DEGs enriched in vector aphids was higher than that in non-vector aphids ( Fig. S2-S5). The COG enrichment analysis showed a similar pattern in that the number of categories in the vector aphids was greater than that in the non-vector aphids ( Fig. S6-S9).

KEGG pathways analysis
To further understand the metabolic pathways that contained differentially expressed genes, KEGG signi cant enrichment analysis was performed on the differentially expressed genes between virustreated vector and non-vector aphids. Among all of the DEGs, there were 287, 74, 91 and 27 DEGs that were assigned to KEGG pathways in Sg-BYDV, Sg-WDV, Rp-BYDV and Rp-WDV, respectively (Fig. S10-S13). The DEGs of the richest and shared pathways were associated with metabolism, including translation and fatty acid and carboxylic acid metabolic pathways. Genes involved in energy production pathways may be important for immune defense responses [27]. In the KEGG pathway analyses, we also identi ed DEGs involved in immune pathways.

Immune-related DEGs of vector and non-vector aphids feeding on virus-infected wheat
Feeding on virus-infected plants resulted in DEGs related to the innate immunity pathway being upregulated. The DEGs related to immunity included several members of the pathway of the lysosomes, the JAK-STAT signaling pathway, antigen processing and presentation, ubiquitin mediated proteolysis, and the peroxisome (Table 4). For the DEG analyses, we identi ed 12 DEGs involved in immune pathways in the vector aphids (Sg-BYDV), the majority of which were upregulated (11 DEGs). For the non-vector aphids, there were fewer DEGs related to immunity, and only one DEG involved in immunity was identi ed in each of the Sg-WDV, Rp-BYDV, and Rp-WDV conditions (Table 4).
Cloning and characterization of STAT5B in R. padi For the immunity related genes analysis, STAT5B (signal transducer and activator of transcription 5B) was the only related gene that was upregulated in R. padi fed on WDV infected wheat. STAT5B is a key gene involved in the JAK-STAT pathway, which is important for the innate immune response of insects. Based on the sequence of the transcriptome, we cloned STAT5B from R. padi (designated as RpSTAT5B). The cloned product was validated using RT-PCR and sequencing. The sequencing results of the RpSTAT5B individual isolated colonies con rmed a single transcript identical to that found in the transcriptome sequencing. We have submitted the sequence of STAT5B to the NCBI GenBank database Multiple protein sequence alignment of RpSTAT5B with six STAT5B proteins from other insects was conducted. Structural analysis showed that the RpSTAT5B protein has the typical structural features of the STAT5B family (Fig. 5). Phylogenetic analysis indicated that the amino acid sequence of the RpSTAT5B protein showed high identity to that of the corresponding proteins from other aphids of Homoptera, particularly Rhopalosiphum maidis (Fig. 6). Both the RT-qPCR and transcriptome results showed that RpSTAT5B was upregulated in R. padi fed on WDV infected wheat (Fig. 4), suggesting that RpSTAT5B is critical to the response to feeding on infected wheat plants. However, the mechanism of RpSTAT5B in aphid immune defense remains to be further studied.

Discussion
Our study investigated the effects of plant viruses on vector and non-vector insects at the transcriptional level. The number of DEGs was greater in the vector aphids (Sg-BYDV with 1525) than in the non-vector aphids (Sg-WDV, Rp-BYDV, and Rp-WDV with 494, 589 and 343, respectively). For all DEGs, the analysis showed that many of the DEGs were involved in development, growth and reproduction. The numbers of genes related to these processes were higher in the vectors (Sg-BYDV) than in the non-vectors (Sg-WDV, Rp-BYDV, and Rp-WDV).
For the GO analyses of all DEGs, the genes associated with the cuticles were mostly downregulated.
Particularly in the vector aphids (Sg-BYDV), there were 16 genes associated with cuticles that were differentially expressed, among which 15 genes were downregulated. Moreover, fewer related DEGs were identi ed in the non-vectors, i.e., two DEGs in Sg-WDV, ve DEGs in Rp-BYDV, and two DEGs in Rp-WDV. The cuticle is of great importance for the survival of insects, since it is the articulated exoskeleton that protects the body against the invasion of pathogens [28,29]. The cuticle and its primary biopolymer components, cuticular proteins and chitin, of insects are periodically turned-over and new cuticle is secreted by the insect epidermis during ecdysis (molting) to accommodate the rapid growth and expansion of the body, and it is thought that temporally and spatially-dynamic epidermal expression of diverse cuticular and endocuticle proteins occurs to support the structure of different hard and soft cuticles of insect body parts during development [30]. The down-regulation of genes associated with the cuticles in aphids under wheat virus infection may signify a delay in molting or turnover of cuticleassociated proteins. One study reported a similar downregulation of six of seven cuticular proteins(CP) transcript sequences of Plodia interpunctella (Indian meal moth) 24 h after exposure to the baculovirus Plodia interpunctella Granulosis Virus (PiGV) [31]. The authors hypothesized that infection suppressed activities of cuticular proteins embedded in the peritrophic matrix, a structural barrier to pathogen attack [32,33]. However, they alternatively offered the possibility that the pathogenic virus may also negatively affect molting.
Insects rely on their immune system to ght against pathogens [34]. After feeding on virus-infected wheat plants, the DEGs related to immunity in S. graminum and R. padi were upregulated, including the JAK-STAT signaling pathway, the lysosome, antigen processing and presentation, ubiquitin mediated proteolysis and the peroxisome. As shown in our results, the number of DEGs related to immunity in S. graminum exposed to BYDV was higher than that in S. graminum exposed to WDV and in R. padi exposed to BYDV and WDV. These results suggest that feeding on virus-infected plants has a greater effect on the immune system of vector insects than that of non-vector insects. The DEGs involved in the cytoskeleton were also differentially expressed, which may be related to the immune response [35]. There have been previous studies showing that viruses can interact with and reorganize host cytoskeleton components for intercellular tra cking and infection processes [36][37][38][39]. In addition, the cytoskeleton is also commonly involved in the intracellular transport of viruses [40][41][42][43].
Lysosome pathway is important in humoral immune response. We identi ed CTSB (cathepsin B), LGMN (legumain), LAMAN (lysosomal alpha-mannosidase), NPC2 (Niemann-Pick C2 protein) associated with lysosome were up-regulated in S. graminum exposed to BYDV; LAMAN was up-regulated in S. graminum and R. padi exposed to WDV. Cathepsins are proteases involved in protein degradation, apoptosis, and signaling, and they regulate viral infection and transmission [48][49][50]. In the green peach aphid, a lysosomal cathepsin B is upregulated following acquisition of the circulative-transmitted virus, Potato leafroll virus (PLRV), and together, the protein and virus colocalize at the cell membranes of midgut cells [50].
The JAK-STAT signaling pathway is important for innate immunity and antiviral responses in insects [51][52][53][54]. This pathway consists of the ligand unpaired (UPD), the receptor Domeless, JAK, STAT, protein inhibitor of activated STAT (PIAS), suppressor of cytokine signaling (SOCS), and signal transducing adaptor molecule (STAM) [44,55,56]. We found that PIAS1 (Sg55104) was up-regulated in the Sg-BYDV aphids; this gene encodes E3 SUMO-protein ligase PIAS1. A previous study has shown that PIAS1 is a speci c inhibitor of Stat1-mediated gene activation [57]. In addition, STAT was up-regulated in the Rp-WDV aphids. JAK activation occurs upon ligand-mediated receptor multimerization, and the activated JAKs subsequently phosphorylate additional targets, including both the receptors and the major substrates, STATs. Phosphorylated STATs enter the nucleus by a mechanism that is dependent on importin a-5 (also called nucleoprotein interactor 1) and the Ran nuclear import pathway. Once in the nucleus, dimerized STATs bind speci c regulatory sequences to activate or repress transcription of target genes (Fig. S14) [58].
STAT5B was the only related gene that was differentially expressed in R. padi fed on WDV infected wheat. RpSTAT5B possesses the SH2 (Src-homology Domain), STAT5_CCD (Coiled-coiled Domain), and the STAT_DBD (STAT_bind) conserved domains, which are the characteristic conserved domains of STAT5 (Fig. 5) [59][60][61][62][63]. For prediction of its three-dimensional structure, the SWISS-MODEL template library (SMTL version 2019-06-27, PDB release 2019-06-21) was searched with BLAST and HHBlits for evolutionary related structures matching RpSTAT5B. The predicted three-dimensional structure of RpSTAT5B is shown in Fig. 7, and its identity with the template (SMTL ID: 1y1u.1.A) is 45.47%. Overall, the identi cation of RpSTAT5B, and its increased expression in R. padi fed on WDV infected wheat, suggests it may ful ll a speci c physiological role. Future studies are needed to understand the function of RpSTAT5B in R. padi fed on WDV infected wheat.
Insect pests and viral diseases are important factors affecting wheat yield. In this study, we analyzed the response of vector/non-vector aphids to wheat virus at the transcriptional level. Among all DEGs analyzed, we identi ed several genes associated with immunity, growth, development, reproduction and the cuticle of aphids. The results revealed that the number of DEGs and the functions related to the DEGs were both higher in vector aphids than that in non-vector aphids. This may imply that wheat viruses have a larger impact on vector than non-vector aphids. Exposure to the virus activated the immunity response of the aphids, particularly vector aphids. These results will provide a reference for investigation of new methods to improve the e ciency of pest control, including prokaryotic expression and gene silencing.

Insect rearing and plant infection by the viruses
Aphids from each of the species S. graminum and R. padi was collected from Yangpingguan (34.85°N, 105.63°E), Shaanxi Province, and they were maintained on wheat seedlings in growth chambers kept at 20 ± 1°C, 65% ± 5% relative humidity and a photoperiod of 16 : 8 h (L : D). Aphids used in the study were parthenogenetic descendants from a single isolated female. BYDV (isolate BYDV-GAV) and WDV were collected from Yangling (34.28ºN, 108.22ºE) and Hancheng (35.47°N, 110.45°E), Shaanxi Province, China, from winter wheat plants. The viruses were identi ed using RT-PCR and agarose gel electrophoresis using the methods brie y described below.
We isolated total RNA from the wheat using TRIzol reagent (TaKaRa, Japan) following the manufacturers' instructions to test whether wheat was infected by BYDV. The cDNA was synthesized using PrimeScript RT Master Mix (TaKaRa, Japan). We isolated total gDNA from the wheat using DNAsecure Plant Kit (TIANGEN, China) following the m anufacturers' instructions to test whether wheat was infected by WDV. Then PCR was performed; the primers used for the virus identi cation are shown in Table S1. The PCRs consisted of the following: 6.5 μl of Reaction Mix (TIANGEN, China), 1 μl of each forward and reverse primer, 1 μl of cDNA/gDNA and 4.5 μl of ddH2O. The cycling conditions were 94°C for 3 min; 35 cycles of 94°C for 30 s, 55°C for 30 s, and 72°C for 1 min; and nally, 72°C for 5 min. Then, agarose gel electrophoresis was performed, and we sequenced the PCR products to con rm that the wheat was infected by BYDV or WDV. The detailed methods of the BYDV-infected and WDV-infected plant treatments, as transmitted by the leafhopper and S. graminum, respectively, have been previously described as follows [64]. Wheat plants at the second leaf stage were rst covered with transparent, plastic, tube-shaped cages (30 cm in height, 13.5 cm in diameter, and with a mesh screen cover on the top), and then 5 late instar nymphs of the vector were introduced into the cage and allowed to feed on the leaves. All the nymphs were removed after 5 days, and the treated plants were maintained under the growth chamber described as above.

Feeding assays and RNA isolation
The young rst instar nymphs of S. graminum and R. padi were fed on BYDV-infected, WDV-infected, and healthy wheat plants to achieve six experimental treatments: the vector treatment was de ned as S. graminum fed on BYDV-infected wheat (Sg-BYDV); the non-vector treatment was de ned as S. graminum fed on WDV-infected wheat and R. padi fed on BYDV-GAV-infected, WDV-infected wheat (Sg-WDV, Rp-BYDV, Rp-WDV); and the control treatments were S. graminum and R. padi fed on healthy wheat seedlings (Sg-ck, Rp-ck).
Three biological replications of each treatment were conducted and each replicate contained 20 individual aphids. For each treatment-replicate, they were maintained in a growth chamber under the same environmental conditions as described above. We collected 20 adult aphids into a 1.5 ml tube, and they were immediately ash frozen in liquid nitrogen and stored at -80 °C. Total RNA was extracted from each sample using TRIzol reagent (TaKaRa, Japan) following the manufacturers' instructions. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using a Qubit® RNA Assay Kit in a Qubit®2.0 Fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Transcriptome sequencing and analysis
A total amount of 3 μg of RNA per sample was used as input material for the RNA sample preparations.
Sequencing libraries were generated using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer's recommendations. Brie y, mRNA was puri ed from the total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out, and the rst and second strands were synthesized. The puri ed and adaptor-ligated cDNA was subjected to PCR ampli cation. Finally, the PCR products were puri ed (AMPure XP system), and the library quality was assessed on the Agilent Bioanalyzer 2100 system.
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2000 platform and paired-end reads were generated. Clean data were obtained by removing reads containing adapters, reads containing poly-N and low-quality reads. At the same time, the Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All of the downstream analyses were based on clean data with high quality. Transcriptome assembly was accomplished using Trinity (r20131110) with min_kmer_cov set to 2 by default and all other parameters set to default [65].
Differential expression analysis of the two groups was performed using the DESeq R package (1.10.1).
DESeq provides statistical routines for determining differential expression of digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted Pvalue <0.05 found by DESeq were considered differentially expressed.
Gene ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the topGO R packages based on the Kolmogorov-Smirnov test. We used KOBAS software to test the statistical enrichment of the differentially expressed genes in the KEGG pathways. Genes with similar expression patterns are usually functionally related. Cluster 3.0 software was used to analyze differentially expressed genes by hierarchical clustering with Euclidean distance as the distance matrix. The clustering results were displayed by Java Treeview.

RT-qPCR validation
To validate differential expression in response to the viruses, 13 transcript sequences were selected in S. graminum and R. padi and then compared with the expression levels between the virus-exposed and nonvirus-exposed aphids. Actin was used as an internal reference gene to normalize the expression level. The primers were designed using Primer Premier 5.0 software and are shown in Table S2. The total RNA extraction and cDNA synthesis methods were performed as described above. Total RNA was isolated from three biological replicates, and each targeted gene included three technical replicates.
The reactions of real-time RT-qPCR were performed in 25 μl volumes using TB Green Premix Ex Taq (TaKaRa, Japan) with 12.5 μl of TB Green Premix Ex Taq II, 1 μl of each forward and reverse primer (10 μM), 0.5 μl of ROX Reference Dye II, 2 μl of cDNA and 8 μl of nuclease free water. The cycling parameters were as follows: 95°C for 30 s and 40 cycles of 95°C for 5 s and 60°C for 30 s. A standard curve experiment was set up by performing a dilution series with 5 -1 dilution times and 5 dilution points to evaluate the PCR e ciency (E). The Q-PCR was performed on a QuantStudio 3 Real-Time PCR System (Applied Biosystems), and QuantStudio™ Design and Analysis Software was used to analyze the PCR assays. The relative quantitative data analysis was performed using the 2 -△△ct method.
Cloning and characterization of the STAT5B STAT5B (signal transducer and activator of transcription 5B) was the only relevant gene that was differentially expressed in R. padi fed on WDV infected wheat. To understand the mechanism of STAT5B in the insect immune defense response, we designed primers using Primer Premier 5.0 to clone STAT5B based on the sequence found in the transcriptome (Table S3). PCR analyses were conducted on a PCR Ampli er (purchased from Eppendorf, Yangling, China) using Tks G ex™ DNA Polymerase (purchased from TaKaRa, Yangling, China). Cloned products were validated by RT-PCR and sequencing. Multiple sequence alignments were performed using Clustal W and DNAMAN. The phylogenetic tree was generated in Mega 6 using the Neighbor-Joining method. All sequences used in the analysis are listed in Table S4. The protein structure model was predicted using WISS-MODEL (https://swissmodel.expasy.org/) and visualized in Pymol (v1.3r1).