Author Correction: Transcriptome analysis of the salivary glands of the grain aphid, Sitobion avenae

An amendment to this paper has been published and can be accessed via a link at the top of the paper.

The salivary components of insects are thought to play crucial roles during interactions with host plants. During feeding, insects secrete saliva with wide composition into host plants for ingesting nutrients and degrading toxins 1,2 . Saliva also contains some proteins and small molecules that have been discovered to activate or suppress plant defense responses [3][4][5][6] .
Unlike leaf chewing insects such as lepidopterans, aphids (Hemiptera: Aphidoidea) have highly modified piercing-sucking mouthparts (stylets) that can penetrate between plant cells to feed phloem sap from sieve elements 7 . During probing and feeding, aphids secrete two types of saliva: gelling saliva, which solidifies into a tube-like sheath to protect the stylet from mechanical damage and chemical attack, and watery saliva, which is secreted into plant cells, intercellular spaces and phloem [8][9][10] . The feeding process of aphids is similar to the infection process of plant pathogens, with an interplay between aphid and host plants that follows the plant-pathogen model proposed by Jones and Dangl 11,12 . Pathogens deliver effectors into the host to modulate plant immunity using a specific secretion system 13,14 . Aphids are also thought to secrete some salivary proteins as effectors into their host plants to alter cell metabolic processes and modulate plant defense responses 15,16 .
Recently, a number of studies have focused on the function of salivary proteins during aphid-host interactions. Salivary proteins between 3 and 10 kD of the green peach aphid, Myzus persicae, can elicit plant defense responses in Arabidopsis thaliana 17 . Salivary secretions by aphids can also prevent phloem from clogging as a result of a plant wound response to enable phloem feeding to continue 18 . The essential roles of salivary protein C002 in the successful feeding of the pea aphid, Acyrthosiphone pisum, on fava bean have been demonstrated using RNA interference 19 . Overexpression of M.persicae effector MpC002, potato aphid Macrosiphum euphorbiae effectors Me10 and Me23 in Nicotiana benthamiana increased aphid fecundity, suggesting their ability to suppress plant defenses, however, in planta expression of Mp10 and Mp42 reduced aphid fecundity and Mp10 induced obvious chlorosis response in N. benthamiana, suggested these effectors involved in the activation of plant defense response [20][21][22] .
For investigating the roles of saliva in aphid-plant interactions, the composition of aphid salivary proteins first needs to be identified. Actually, the transcriptome of salivary glands has been analyzed in some hemipteran species such as potato leafhopper Empoasca fabae 23 , whitefly Bemisia tabaci 24 , rice brown planthopper Nilaparvata lugens 25 , green rice leafhopper Nephotettix cincticeps 26 , tarnished plant bug Lygus lineolaris 27 . Some salivary proteins with enzymatic activities were also identified in the watery saliva of M. persicae using mass spectrometry (LC-MS/MS); e.g., glucose oxidase, first identified in the chewing insect Helicoverpa zea, was demonstrated to be an effector that suppresses plant defense 4,28 . In a dual transcriptomic-proteomic approach, over 300 secretory salivary proteins from the salivary glands of A. pisum were identified 29 .
The grain aphid, Sitobion avenae, is an important agricultural pest of cereals causing serious economic losses through nutrient robbing and transmitting plant viruses (BYDV) 30,31 . In research to elucidate the composition and function of S. avenae saliva, pectinase was detected, and exogenous application of pectinase induced volatile emissions in wheat and attracted the aphid parasitoid Aphidius avenae 32 . Ma et al. found polyphenol oxidase (PPO) in S. avenae saliva, and wheat seedlings treated with PPO had increased expression of genes related to plant defense signaling 33 . Recently, the composition of watery salivary proteins of S. avenae was studied by tandem mass spectrometry, and 12 proteins were identified 34 .
Aphid salivary secretory proteins were usually collected from artificial diet for further identification by proteomics, some proteins with low concentrations or only induced during interacted with host plants may not be detected. Aphid salivary glands are paired and composed of two primary glands and two accessory glands, secreted salivary proteins are mainly synthesized in the primary glands 35 . Although the function of the accessory glands has not been well studied, they might be involved in the transmission of viruses such as BYDV in S. avenae 36 . Therefore, in this study, we directly dissected whole salivary glands from aphids ( Fig. 1) and sequenced the transcriptome of the salivary glands of S. avenae. Of 33,079 unigenes identified in the whole salivary glands, 526 unigenes were predicted to encode secretory proteins, and some of their orthologs have been proved to play important roles in aphid-host interactions. These findings firstly provide insight into the identification of potential effector molecules in S. avenae saliva and further our understanding of the roles of saliva in aphid-wheat interactions.

Results and Discussions
Illumina sequencing and unigene assembly. In total, 58,224,912 bp and 50,854,016 bp raw reads were acquired from salivary glands of apterous and alate S. avenae adult, respectively. After removing adapters, ambiguous nucleotides and low quality sequences, 56,280,842 bp and 49,104,484 bp cleans reads remained. Subsequently, the transcriptome of S. avenae salivary glands was de novo assembled using the short reads assembling program-Trinity 37 , which were then clustered into 41,335 transcripts and 33,079 unigenes (Table 1). These transcripts ranged from 201 to over 14,923 bp with an average size of 941 bp. Among the transcripts, 22,532 (54.51%) were between 200 bp and 500 bp long, and 5,510 (13.33%) were over 2,000 bp. Also, among the assembled unigenes, 20,698 (62.57%) unigenes were between 200 bp and 500 bp long, and 3,120 (9.43%) were over 2,000 bp; mean length of unigenes was 711 bp (Fig. 2). All sequences of the unigenes in this study are provided in (Supplementary Table 1
EuKaryotic Orthologous Groups (KOG) is a version of the Clusters of Orthologous Groups (COG) for identifying orthologous and paralogous proteins in eukaryotic organisms 39 . In total, 7,956 unigenes were categorized into 26 groups, among these categories, the cluster of "General function prediction only" group had the most unigenes (1,425, 17.89%) followed by "Signal transduction mechanisms" (1,186, 14.91%), "Posttranslational     Figure 1).

Metabolic pathway analysis by Kyoto Encyclopedia of Genes and Genomes (KEGG). When the
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was used to describe the network of molecular interactions and metabolic pathways in cells 40 , 6,772 unigenes in the salivary glands of S. avenea were mapped to a total of 229 KEGG pathways. Among these pathways, "Ribosome" (384 unigenes), "Protein processing in endoplasmic reticulum" (175 unigenes), and "PI3K-Akt signaling pathway" (166 unigenes) had the most unigenes (Fig. 4).
Putative secretory proteins. Salivary proteins are proposed to be involved in interactions with plant, only when they can be secreted into plant during aphid probing and feeding. Therefore, all unigenes were analyzed for the presence of signal peptide and potential cleavage site using SignaIP software in this study. In total, 526 putative secretory proteins were obtained from all unigenes, however, some putative secretory proteins might have been missed as a result of partial ORF sequences, such as 3′ partial, 5′ partial and internal partial sequences. Of all putative secretory proteins, 335 (63.69%) proteins were functionally annotated and 191(36.31%) showed no similarities with known-function proteins in the Nr database. And, some function-annotated putative secretory proteins are very closely related to several insect salivary proteins that have been proved to play important roles in insect-host interactions, such as digestive and detoxifying enzymes and effectors (Supplementary Table 2).
Digestive and detoxifying enzymes. During phloem feeding, aphids secrete salivary proteins with digestive enzyme activity that facilitate probing and feeding. The transcript of a cell wall degradation enzyme, beta-mannosidase (c15158_g1) was found in our study. Beta-mannosidase catalyzes the endo-wise hydrolysis of the backbone of mannan and heteromannans, including a major component of plant cell walls, hemicellulose polysaccharides 41 . This enzyme can help aphids penetrate the plant cells. Sugar-degrading enzymes and proteases such as maltase, beta-glucuronidase, serine protease, trypsin and cathepsin were also identified in the salivary glands of S. avenae in our study. Wheat phloem sap is high in sucrose and some amino acids, predominantly glutamic acid, aspartic acid and serine 42 . The presence of these secretory digestive enzymes from S. avenae could thus aid extra-orally digestion, consistent with the feeding style of S. avenae.
Peroxidase (c28134_g1, c13443_g1, c12708_g1, c14634_g2, c14243_g1, c14995_g1), an oxidoreductase was detected here, and was reported previously in salivary glands of Hession fly Mayetiola destructor 53 , A. pisum 29 and M. viciae 51 . Reactive oxygen species (ROS), including singlet oxygen, superoxide and hydrogen peroxide are involved in signaling pathways to activate plant defense responses and resistance to aphid. For example, accumulation of hydrogen peroxide at aphid feeding sites indicated reactive oxygen species are involved in early signaling of A. thaliana after infestation by the cabbage aphid Brevicoryne brassicae 54 . High levels of H 2 O 2 accumulated in a resistant near-isogenic wheat line but not in a susceptible line after infestation by the D. noxia 55 . A sudden rise in oxygen concentration in sieve elements might result in protein coagulation in cucurbit phloem sap 56 . Additionally, high H 2 O 2 levels in plants could be toxic to aphids. Peroxidase, acting as antioxidant enzymes involved in H 2 O 2 scavenge, suggest that peroxidase detected in aphid salivary glands may protect aphids from plants oxidative stress as a detoxifying enzymes and play important roles in suppressing ROS production and ROS-induced plant defense responses.
Effectors eliciting or suppressing plant defenses. We found beta-glucosidase (c10709_g1) in the salivary glands of S. avenae. Beta-glucosidase, mainly regarded as a digestive enzyme for molecules such as hemicellulose, was also detected in the salivary glands of the termite Neotermes koshunensis 57 . However, it might also be an effector to activate plant defense responses. Lima bean Phaseolus lunatus treated with a solution of beta-glucosidase emitted more volatile compounds, which were similar to those emitted in response to red spotted spider mite Tetranychus urticae infestation 58 . Leaves treated with commercial beta-glucosidase released volatile blends similar to that of leaves infested with cabbage white butterfly Pieris brassicae 59 . Also, the levels of plant defense signaling molecules salicylic acid, ethylene, and H 2 O 2 in rice increased after the application of beta-glucosidase, and the defense signaling pathways induced by beta-glucosidase were similar to those activated by infestation with N. lugens 60 .
Lipases (c11727_g1, c8374_g1) and some phospholipases (c13818_g1, c7165_g1) were detected in the transcriptome of S. avenae salivary glands in this study. Lipases also have been found in the salivary glands of some insects, including mosquito Anopheles stephensi 61 , M. destructor 62 , large milkweed bug Oncopeltus fasciatus 63 and E. fabae 23 . The primary function of lipases is presumed to be the breakdown membrane lipids and thereby the membrane; thus, they were first predicted to act as digestive enzymes to facilitate penetration of the cell membrane. However, Schäfer et al. showed that lipases in grasshopper oral secretions induced accumulation of cyclopentenone 12-oxo-phytodienoic acid (OPDA), a precursor of JA biosynthesis in A. thaliana, and external application of lipase solution to wounded leaves also highly increased the levels of 13-hydroperoxy SCIEntIFIC RepoRTS | 7: 15911 | DOI:10.1038/s41598-017-16092-z octadecatrienoic acid, OPDA, JA, and jasmonic acid-isoleucine, suggesting that lipases elicited plant defense responses 64 . Phospholipases hydrolyze phospholipids, which take part in lipid synthesis, lipid-derived signaling pathway and plant stress responses. It is assumed that phospholipases D (PLD) and its products, phosphatidic acid, are involved in plant signal transduction cascades and the lipid metabolic pathway to influence plant stress responses. Gene expression and enzyme activity of PLD increased rapidly after exposure to various stresses, such as mechanical wounding, frost and pathogen infestation [65][66][67] . Increases in the PLD transcripts and accumulation of PLD along the plasma membrane were observed after rice leaves were infected with the bacterial pathogen Xanthomonas oryzae pv. oryzae 68 . Also, PLD participates in the production of plant defense response signaling molecules (ABA, ET, NADPH oxidase) 66 . Thus, lipases in S. avenae or lipase-derived molecules are likely to serve as effectors to induce downstream plant defense responses.
Salivary glucose oxidase (GOX) was the first effector identified in the saliva of herbivores to suppress plant defense. The caterpillar H. zea secretes GOX into its host Nicotiana tabacum to suppress nicotine production and the jasmonic acid defense signaling pathway 4 . GOX was found in the saliva of many other caterpillar species and in the aphid M. persicae 28 , but not in other aphid species. We identified another highly expressed glucose-methanol-choline oxidoreductase, glucose dehydrogenase (GLD, c10172_g2, c10815_g1, c12301_g1), in the salivary glands of S. avenae. GLD is also found in the saliva and salivary glands of other aphid species 29,34,50,51 . We predicted that GLD functioned in a similar way to GOX and as a potential effector to suppress plant defense responses after aphid infestation.
C002 is an aphid-specific watery saliva protein and also a well-known effector protein related to aphid feeding behavior and subsequent survival and fecundity [69][70][71] . Knockdown of C002 transcript of A. pisum and S. graminum resulted in high mortality of aphids 19,69,71 . The reproduction rate of M. persicae increased after feeding on host plants that over-express MpC002 but decreased after feeding on plants producing double-strand RNA (dsRNA) against C002 20,70 . At the proteomic level, C002 protein has been identified in M. persicae and A. pisum saliva 19,28,72 . However, its underlying mechanism is still unknown. An ortholog of C002 (c12732_g1) was also found in the salivary glands of S. avenae in this study with very high abundance in the transcriptome (RPKM = 10,926.33).
Angiotensin converting enzyme-1 (ACE-1, c10308_g1) was found in the salivary glands of S. avenae in the present study. ACE is a zinc-metallopeptidase found on the endothelial, epithelial and neuronal tissues in mammals 73 . The function of ACEs were studied in various insects such as locust Locusta migratoria 74 , cotton leafworm Spodoptera littoralis 75 and buffalo fly Haematobia irritans exigua 76 . Three ACE genes, ACE1, ACE2 and ACE3, were identified in the genome of A. pisum 77 , and the ACE proteins were also detected in the saliva of A. pisum 29,72 . ACE1 and ACE2 are highly expressed in the salivary glands, and compared with aphid survival rate on artificial diets, the survival rate of aphids with simultaneous knockdown of ACE1 and ACE2 decreased significantly after aphids fed on plants. These results indicated that ACE1 and ACE2 can function as effectors to modulate plant physiological processes to benefit aphid infestation 77 .
Trehalose, a kind of glucosidase in many organisms ranging from bacteria, fungi to plants 78 . Although the level of trehalose is very low in plants, many studies have demonstrated that it plays a regulatory role in sugar metabolism, growth, development and stress responses of plants 79 . Overexpressing trehalose biosynthetic genes of microbial and plant origin or receive exogenous trehalose, stress tolerance increased in several plants, such as tobacco, potato and rice 80,81 . We identified trehalase (c11521_g1, c13987_g1, c14644_g1), which degrades trehalose, in the salivary glands of S. avenae, and it has also been found in salivary glands of B. tabaci and the saliva of M. dirhodum and S. avenae 24,34 . Trehalase activity in A.thaliana increased before trehalose accumulated when infested with the trehalose-producing pathogen Plasmodiophora brassicae, suggesting that trehalase may be a part of the plant defense responses and prevent excess accumulation of trehalose in plant cells 82 . Aphid-secreted trehalase may also play a role in interrupting trehalose accumulation to repress plant defense responses. However, a study showed that exogenous trehalose suppressed the transcript levels of some genes encoding enzymes related to plant defense, wound response, or pathogenesis, such as peroxidase-2 (PRXR2), basic endochitinase (ChiB), endo-1,3-b-d-glucanase (BGL1), lipoxygenase-2 (LOX2), and a chitinase-like protein 1 (CTL1) 83 . The role of trehalase in plant defense response needs further study.
Among the putative secretory proteins, we found odorant binding proteins (OBPs, c5921_g1, c2376_g1) and chemosensory proteins (CSPs, c5314_g1, c8576_g1, c18419_g1, c2724_g1). Transcripts of OBPs and CSPs have also been found in salivary glands of N. lugens, M.persicae, M.cerasi and bird cherry-oat aphid Rhopalosiphum padi 25,84 . Insect OBPs and CSPs involved in olfaction and gustation are thought to be crucial for insect behaviors such as locating food and ovipository sites, as well as intraspecific communication 85 , and are mainly specifically expressed in chemosensory organs such as antennae and mouthparts and predicted to function in chemoperception [86][87][88] . However, some insect OBPs and CSPs are found in other tissues such as legs, heads, bodies and salivary glands 86 with functions in insect development 89 , leg regeneration 90 , immune responses, and even interactions with the host. For example, Dengue virus (DENV) infection increased transcripts expression of OBP10 and OBP22 in the salivary glands of Aedes aegypti and silencing of the OBP10 or OBP22 genes resulted in a low efficiency of mosquito blood-feeding 91 . OBPs were found in the salivary glands of Anopheles gambiae and predicted to be secreted into host cells to manipulate host physiology by scavenging host amines 92 . Also, CSPs were identified using mass spectrometry as the most abundant proteins in the mandibular glands of larvae of the butterfly Vanessa gonerilla and speculated to play important roles in detecting microorganisms on plant surfaces, recognizing host plants and communicating with conspecifics 93 . CSP4 in cotton bollworm Helicoverpa armigera and oriental tobacco budworm Helicoverpa assulta possess unique functions that act as surfactants to reduce water surface tension and, consequently, pressure during sucking 94 , suggesting CSPs are involved in insect feeding. Among the CSPs identified in S. avenae, unigene c18419_g1 had high similarity to Mp10 known to be an effector in M. persicae. Mp10, a chemosensory protein was detected in the heads, salivary glands and whole bodies of M. persicae. Functional assays showed that overexpression of Mp10 in N. benthamiana suppressed bacterial pathogen-associated molecular pattern (PAMP) flg22-induced defense responses, but induced chlorosis and local cell death in N. benthamiana, resulting in a decrease of M. persicae fecundity 20,22 . Thus, OBPs and CSPs secreted by salivary glands of S. avenae may play important roles in aphid-host interactions and warrant further functional study.
Among potential effectors that were detected, such as lipid-binding proteins apolipophorins (c26673_g1) 51 , we also identified transcript sequence c10120_g1 as an ortholog of potato aphid M. euphorbiae effector Me10. Overexpression of Me10 in the host plant N. benthamiana increased M. persicae fecundity, suggesting its ability to suppress plant defenses 21,52 . As we mentioned before, some putative secretory proteins with effector activity may be missed because of partial sequences. For example, the unigene c9478_g1 (5′ primer partial) is a likely homolog of a potential salivary effector Mp55 of M. persicae, A. thaliana expressing Mp55 increased aphid reproduction, and accumulated less 4-methoxyindol-3-ylmethylglucosinolate, callose and hydrogen peroxide in response to aphid infestation, suggesting a role of Mp55 in suppressing plant defenses 95 .
Calcium ions (Ca 2+ ) constitute a ubiquitous intracellular second messenger in many plant signaling pathways including induction of defense responses 96,97 . Mechanical damage of phloem tubes can trigger the sieve plate occlusion to avoid the outflow of phloem sap because of the releasing of Ca 2+ into the sieve element lumen 98 . But damage caused by aphid stylets penetration don't lead to the phloem sieve cell plugging 18 . Will T. firstly reported that M. viciae saliva contained some Ca 2+ binding proteins prevent Ca 2+ -dependent sieve occlusion in Vicia faba through inducing dispersed forisomes, a Ca 2+ -driven contractile protein that can cause reversible plugs in sieve element, to return to the non-plugging state, and secretion of watery saliva also seems to be a universal way for aphids to suppress sieve-plate occlusion 99,100 . Some Ca 2+ binding proteins were also found in salivary glands and saliva of hemipterans using transcriptomic and proteomic technologies, such as regucalcin in A. pisum and N. cincticeps 26,72 , the most highly expressed salivary glands gene calcium-binding protein SP84(NcSP84) in N. cincticeps 26 . These indicated that Ca 2+ binding proteins secretion to suppress plant defense may be a common strategy amongst phloem-feeding insects.
Interestingly, only 4 transcripts were annotated in top 15 highly expressed unigenes, with most still categorized as function unknown. Among these 4 annotated unigenes, GLD (c12301_g1, c12029_g2) was described before and predicted to be an aphid effectors to promote aphid infestation. Unigene c13074_g1 and c29740_g1 were annotated as MLH and TPR, respectively. MLH is a DNA binding protein and involved in plant salt tolerance mechanisms. The gene expression for MLH polyproteins is upregulated in the salt-tolerant genotype of Acer palmatum, suggesting it is involved in the improvement of plant resistance to abiotic stress 101 . The TPR motif is a protein-protein interaction module, which is important for the functioning of chaperones, cell-cycle, transcription, protein transport complexes and the gibberellin signal transduction pathway 102, 103 . Gibberellins (GAs) represent an important class of plant hormones that control growth and developmental processes. We predict that TPR secreted by aphids may have roles in the interactions with plant proteins and affect plant development. The RT-PCR and RT-qPCR results showed that all of these highly expressed transcripts were specifically expressed in aphid salivary glands (Fig. 5a,b). Relative expression of the most highly expressed unigenes c10120_g1 and c12029_g2 (GLD) in salivary gland was 238.5 ± 27.7 and 823.3 ± 51.9 times higher than in the whole body, respectively. These highly expressed salivary gland-specific genes may play vital roles in aphid-plant interactions, and are worth further analysis.
Among 9 genes of interest, β-MAN, MAL and PLD were detected in the transcriptome of S. avenae salivary glands, but the RT-PCR and RT-qPCR results showed that the expression levels of these transcripts in alimentary canal were significantly higher than in other tissues. The expression levels of ACE and putative Ca 2+ binding proteins RCN in salivary glands were significantly higher than those in other tissues. GST-1 transcript was specifically expressed in the salivary glands. Aphid salivary protein with GST activity involved in modifying plant defense responses, overexpression of a single putative GST protein named Me47 detected in M. euphorbiae saliva could induce or suppress plant defense depending on different host plants 104 . TRE and CALU had high expression in both salivary glands and alimentary canal. In our study, RT-PCR and qPCR results showed that RGN was very low in the salivary glands compared with its levels in the alimentary canal and whole body (Fig. 6a,b). RGN was found exclusively in the anterior fat bodies of the flesh fly Sarcophaga peregrine 105 , whereas the RGN protein was detected in the watery saliva of A.pisum 34 , suggesting that RGN may be mainly produced in the alimentary canal and other tissues such as fat bodies, then transported from the haemolymph into the salivary gland to be secreted into the host.
In conclusion, we revealed the transcripts of the salivary glands of S. avenae using Illumina HiSeq 2500. Five hundred and twenty five putative secretory proteins that were expected to be secreted into plants and may play critical roles in aphid-host interactions were identified. Among of them, some highly and salivary gland-specifically expressed genes were also uncovered, but most of these putative secretory proteins were still function-unknown, which are worthy of further study. We also suggest a model describing the potential roles of some salivary proteins in aphid-host interactions (Fig. 7). Further investigations are needed to confirm these predicted secreted proteins at the proteomic level for a more comprehensive understanding of the composition of

Materials and Methods
Insect rearing. A clone of S. avenae was initially established from a single aphid collected from wheat field in Langfang, Hebei Province, China and has been reared on wheat plants (variety Beijing 837, which is susceptible to S. avenae) for 5 yr in an indoor environment at 20 ± 1 °C, 75-80% relative humidity and 16 h light/8 h dark.
Sample collection and RNA isolation. About 600 pairs of salivary glands were individually dissected from apterous and alate adult aphids in phosphate buffered saline (pH = 7.2, Hyclone, Thermo Scientific, MA, USA) respectively, then quickly transferred to TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) on ice. Total RNA was extracted using TRIzol Reagent following the manufacturer's instructions and stored at −80 °C until used. RNA concentration was measured using Qubit RNA Assay Kit and Qubit 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit and the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Three micrograms total RNA sample with standard quality (1.8 < OD260/280 < 2.1, RIN values > 8.0) was prepared for further sequencing.
Illumina sequencing, assembly, and annotation. NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (New England Biolabs (NEB), Beverly, MA, USA) was used to generate the sequencing libraries following the manufacturers' introductions. Firstly, mRNA was purified from total RNA sample using poly-T oligo-attached magnetic beads. Fragmentation was performed using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First strand cDNA synthesis was carried out using M-MuLV Reverse Transcriptase and random hexamers, then second strand cDNA was synthesized using DNA polymerase I and RNase H. The remaining overhangs were converted into blunt ends using exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptors were ligated to DNA fragments for hybridization. The library fragments were purified with the AMPure XP system to select the cDNA fragments with a length between 150 and 200 bp. Three microlitres USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated DNA at 37 °C for 15 min followed by 5 min at 95 °C, PCR was then performed to amplify cDNA with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. The PCR products were purified with AMPure XP system and quantified using the Agilent Bioanalyzer 2100 system (Agilent Figure 7. Schematic drawing of potential roles of secretory proteins from aphid saliva in aphid-plant interactions. During probing and feeding, the aphid secretes saliva into the cytoplasm of the plant cell. Digestive and detoxification enzymes ( ) in salivary proteins facilitate probing and feeding because these enzymes are involved in the breaking down cell walls and membranes, metabolizing sugars and amino acids, and detoxifying secondary metabolites in plant tissues. Some digestive enzymes and degradative products also can induce plant defense responses as potential effectors. Some effectors ( ) in saliva can be recognized by plant transmembrane pattern recognition receptors (PRRs), resulting in the induction of plant basal defenses (pattern-triggered immunity, PTI). Aphid also secretes other effectors ( ) to suppress plant defenses to promote infestation, resulting in effector-triggered susceptibility (ETS). However, in resistant interactions, an effector ( ) can be specifically recognized by plant resistance proteins (R proteins) according to the gene for gene hypothesis, resulting in a stronger defense responses, effector-triggered immunity (ETI) such as induced cell death. Penetration of plant cell membranes by aphids causes an increase in the cytosolic concentration of Ca 2+ ([Ca 2+ ] cyt ), resulting in sieve pore occlusion and calcium signaling for defense responses such as the reactive oxygen species (ROS) burst. Ca 2+ -binding proteins ( ) injected with saliva can bind Ca 2+ , thereby preventing sieve pore blockage and calcium signaling pathways. Also, the hole in the membrane created by the aphid can be sealed after the secreted gelling saliva solidifies into a sheath, minimizing the influx of Ca 2+ and downstream plant defense responses. However, few effectors have been identified, and details of the mechanisms involved in aphid-plant interaction are still unknown. Abbreviations: SEs: sieve elements MAPK: mitogen-activated protein kinase; MeSA: menthyl salicylate; MeJA: methyl jasmonate.
Technologies, CA, USA). The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, China) according to the manufacturer's instructions.
After cluster generation, the library preparations were sequenced using an Illumina HiSeq. 2500/Miseq platform and paired-end reads (the sequencing strategy was PE150). The clean reads were obtained after adaptor sequences, ambiguous "N" nucleotides (the percentage of "N" > 10%) and low-quality sequences (the ratio of nucleotides with Qphred ≤ 5 was more than 50%) were removed from raw reads. The clean reads were assembled using Trinity r20140413p1 min_kmer_cov:2 and the other default parameters as described for de novo transcriptome assembly without a reference genome to generate transcripts and unigenes 109 . For homology searches and annotation, all unigenes were used in a search of public databases including non-redundant protein (Nr, e-value ≦ 1.0e −5 ), nucleotide sequence (Nt, e-value ≦ 1.0e −5 ), Pfam (e-value ≦ 0.01), euKaryotic Ortholog Groups (KOG)/Clusters of Orthologous Groups of proteins (COG, e-value ≦ 1.0e −3 ) and Swiss-Prot (e-value ≦ 1.0e −5 ). Functional annotation by Gene Ontology (GO) term was analyzed using Blast2GO with a cutoff value of E-value ≦ 1.0E −6 . Pathways were annotated by KEGG Automatic Annotation Server (KAAS) based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) with a cutoff value of E-value ≦ 1.0E −10 .
Putative secretory proteins. Open reading frames (ORFs) within transcript sequences generated by de novo RNA-Seq transcript assembly using Trinity were identified using TransDecoder v3.0.0 (https://github.com/ TransDecoder/TransDecoder/releases). Signal peptides and cleavage sites in amino acid sequences were predicted by the SignaIP 4.1 Server (http://www.cbs.dtu.dk/services/SignalP/). For transmembrane domains prediction, amino acid sequences with a signal peptide were submitted to the TMHMM Server v. 2.0 (http://www.cbs.dtu.dk/ services/TMHMM/). Putative proteins with a signal peptide and 0-1 transmembrane domain (the signal peptide can be a transmembrane domain) were considered to be potential secreted proteins 25 .
RT-PCR and RT-qPCR. Total RNA was extracted from 600 salivary glands, 450 alimentary canals, 10 whole bodies of apterous adult aphids and 10 whole bodies minus salivary glands as mentioned before. The cDNA was synthesized from 1 μg RNA using TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China) for RT-PCR and qRT-PCR. All specific primers for RT-PCR and qRT-PCR were designed with Primer Premier 5.0 (PREMIER Biosoft, CA, USA) and are shown in (Supplementary Table 4). β-Actin and NADH dehydrogenase (NADH) were used as reference genes to normalize target gene expression 110,111 . All PCR products were sequenced (Sunbiotech, Beijing, China).
RT-PCR was conducted in a 20 μL reaction volume containing 10 μL 2 × Taq PCR MasterMix (BioMed, Beijing, China), 2 μL cDNA, 1 μL each forward and reverse primer (10 μM), 6 μL ddH 2 O. The RT-PCR conditions were 3 min at 94 °C; 28 cycles of 30 s at 94 °C, 30 s at 60 °C and 30 s at 72 °C; and final 10 min at 72 °C. PCR products were separated in 1% agarose gels and stained with ethidium bromide (Thermo Scientific, MA, USA). RT-qPCR was performed using an ABI 7500 Real-Time PCR System (Applied Biosystems, CA, USA).cDNA was diluted 10-fold and then used as the template to determine the relative expression of the target gene in a 20 μL reaction volume containing 2 μL cDNA, 0.5 μL each of 10 μmol L −1 forward primer and reverse primer, 10 μL of 2 × SYBR premix Ex Taq (Tli RNaseH Plus, Takara, Dalian, China) and 0.4 μL of 50 × ROX Reference Dye II (Tli RNaseH Plus, Takara, Dalian, China) at the following conditions: 30 s at 95 °C; followed by 40 cycles of 30 s at 95 °C, and 30 s at 60 °C. RT-qPCR reaction for each sample was carried out with 3 biological replicates and 3 technical replicates.
RT-qPCR data analysis. Each treatment was performed in triplicate, and the differential expression was calculated using the 2 −ΔΔCT method 112 . The fold-change in expression of the unigenes in different tissues to the whole body of apterous adult aphid was calculated. Results were expressed as Means ± SE. All data were analyzed using SAS 9.1 software (SAS Institute Inc., NC, USA) and the differences among groups were examined using one-way analysis of variance (ANOVA) test. P values less than 0.05 were considered as statistically significant.