Fast Evolution and Lineage-Specific Gene Family Expansions of Aphid Salivary Effectors Driven by Interactions with Host-Plants

Hélène Boulain; Fabrice Legeai; Endrick Guy; Stéphanie Morlière; Nadine E Douglas; Jonghee Oh; Marimuthu Murugan; Michael Smith; Julie Jaquiéry; Jean Peccoud; Frank F White; James C Carolan; Jean-Christophe Simon; Akiko Sugio

doi:10.1093/gbe/evy097

Boulain, Hélène;Legeai, Fabrice;Guy, Endrick;Morlière, Stéphanie;Douglas, Nadine E;Oh, Jonghee;Murugan, Marimuthu;Smith, Michael;Jaquiéry, Julie;Peccoud, Jean;White, Frank F;Carolan, James C;Simon, Jean-Christophe;Sugio, Akiko

2018-06-18 00:00:00

Effector proteins play crucial roles in plant-parasite interactions by suppressing plant defenses and hijacking plant physio- logical responses to facilitate parasite invasion and propagation. Although effector proteins have been characterized in many microbial plant pathogens, their nature and role in adaptation to host plants are largely unknown in insect herbivores. Aphids rely on salivary effector proteins injected into the host plants to promote phloem sap uptake. Therefore, gaining insight into the repertoire and evolution of aphid effectors is key to unveiling the mechanisms responsible for aphid virulence and host plant specialization. With this aim in mind, we assembled catalogues of putative effectors in the legume specialist aphid, Acyrthosiphon pisum, using transcriptomics and proteomics approaches. We identiﬁed 3,603 candidate effector genes predictedtobe expressed in A. pisum salivary glands (SGs), and 740 of which displayed up-regulated expression in SGs in comparison to the alimentary tract. A search for orthologs in 17 arthropod genomes revealed that SG-up-regulated effector candidates of A. pisum are enriched in aphid-speciﬁc genes and tend to evolve faster compared with the whole gene set. We also found that a large fraction of proteins detected in the A. pisum saliva belonged to three gene families, of which certain members show evidence consistent with positive selection. Overall, this comprehensive analysis suggests that the large repertoire of effector candidates in A. pisum constitutes a source of novelties promoting plant adaptation to legumes. Key words: Acyrthosiphon pisum, salivary proteins, host adaptation, positive selection, pest evolution, plant defenses. Introduction fundamental role in promoting their species richness and di- Insects comprise the most diverse group of metazoans, and versiﬁcation (Wiens et al. 2015). Almost half of the currently evidence indicates that the evolution of herbivory has played a known insect species feed on plants (Wu and Baldwin 2010), The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non- commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com 1554 Genome Biol. Evol. 10(6):1554–1572. doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Evolution and Gene Family Expansions of Aphid Salivary Effectors GBE and herbivorous insect groups exhibit faster rates of diversiﬁ- M. persicae growth (Pitino and Hogenhout 2013). This obser- cation compared with nonherbivorous species (Wiens et al. vation supports specialization of orthologous effectors to dis- 2015). A central hypothesis accounting for higher species rich- tinct plant species during aphid divergence. Some salivary ness in herbivorous insects proposes an evolutionary interac- proteins that are known to contribute to aphid plant exploi- tion between plant defense mechanisms and plant tation are expressed in salivary glands (Wang et al. 2015a, exploitation strategies of insects (Janz 2011). Furthermore, 2015b) and a few showed sites under positive selection (Pitino continuous interactions between host plants and herbivorous and Hogenhout 2013; Thorpe et al. 2016). However, a global insects are predicted to make herbivore generalism difﬁcult and comprehensive evolutionary analysis of aphid salivary and constrain a given insect species to one or a few host genes has yet to be reported. species (Forister et al. 2015). Since plants provide not only A catalogue of putative salivary effectors was created for food resources, but also habitats and mating sites to many A. pisum upon completion of the genome sequence (The herbivorous insects, plant specialization may induce divergent International Aphid Genomics Consortium 2010). Combined selection in insect populations at a range of traits that can lead transcriptomics and proteomics produced a catalog of 324 to reproductive isolation and speciation (Peccoud et al. 2010; secreted proteins (Carolan et al. 2011). A small number of Mullen and Shaw 2014). Attempting to unveil the basic effectors predicted in this catalogue were functionally char- mechanisms of insect herbivory provides opportunities to un- acterized and have been shown to be involved in plant-aphid derstand the evolutionary and mechanistic basis of plant spe- interactions (Mutti et al. 2006, 2008; Wang et al. 2015a, cialization by herbivorous insects, in particular, and the 2015b). Signiﬁcant developments in RNA-seq technology diversiﬁcation of metazoan life, in general. and high-resolution mass-spectrometry (MS) provide new im- Aphids (Insecta: Aphidomorpha) are pests of wild and cul- petus to revise the salivary gene catalogue and to deﬁne a tivated plants that directly reduce plant nutrients by ingesting new and expanded set of candidate salivary effector genes for phloem sap and indirectly cause diseases by transmitting plant further analyses. In addition, the genome sequences of two pathogens (Blackman and Eastop 2000; Harris and specialist aphids, the Russian wheat aphid, Diuraphis noxia Maramorosch 2014). Aphids are also excellent subjects for (Nicholson et al. 2015) and the soybean aphid, Aphis glycines host specialization studies. Their clade is composed of approx- (Wenger et al. 2017), together with that of the generalist imately 5,000 species (Blackman and Eastop 2000), and most green peach aphid, M. persicae,(Mathers et al. 2017)have are considered to be plant specialists (Peccoud et al. 2010). been recently published. These data offer the opportunity to Aphid mouthparts are modiﬁed into a rostrum or beak with better understand the evolutionary dynamics of salivary effec- the mandibles and maxillae forming needle-like stylets. tor candidates, in particular, their suspected role in the adap- Aphids secrete gelling saliva during the early stages of feeding tation of aphid lineages to their host plants (Pitino and to form a feeding sheath surrounding the stylets, and then Hogenhout 2013; Rodriguez et al. 2017). The critical effectors secrete watery saliva into various plant cells (Moreno et al. that drive aphid–plant interactions may display peculiar gene 2011). Saliva contains effectors that modulate physiological expression and evolutionary patterns, which we search by responses to herbivory and permit feeding (Rodriguez et al. combining transcriptomics and proteomics in A. pisum and 2017). These effectors are likely exposed to natural selection, by conducting a comparative analysis of aphid genomes. in particular by plant surveillance systems and defense mech- anisms (Will et al. 2013). Interference with plant defenses Materials and Methods through various mechanisms has been demonstrated for sev- Aphids, Plants, and Growth Conditions eral effectors secreted by microbial plant pathogens, which ultimately promotes persistence and even spread of these Acyrthosiphon pisum lineage LSR1 (used for whole-genome pathogens (Varden et al. 2017). A subset of effectors, the sequencing; The International Aphid Genomics Consortium so-called avirulence proteins, are detected by plant surveil- 2010) was maintained in a growth chamber at 18 Cwith lance systems and trigger strong immunity in speciﬁc plants, a 16 h day/8 h night photoperiod on broad bean, Vicia faba determining the incompatibility (Bent and Mackey 2007). (Castel cultivar), at low density to avoid the production of Effector genes are diverse, making prediction of effector winged individuals. Vicia faba was grown in a growth cham- functions often difﬁcult from the amino acid sequences. As a ber at 18 C with a 16 h day/8 h night photoperiod for 10 days result, relatively few salivary effectors from aphids have char- before installation of the aphids. acterized interactions with active host defense responses. In planta expression of salivary effectors C002, Mp1,and Mp2 RNA Sequencing from the green peach aphid, Myzus persicae, increases M. persicae fecundity on the host plants Arabidopsis thaliana To prepare RNA samples from aphid salivary glands and ali- and Nicotiana benthamiana, whereas expression of ortholo- mentary tracts, 9-day-old individuals reared at a density of gous genes from another aphid species (the pea aphid, 10–15 aphids per V. faba plant were rapidly dissected with Acyrthosiphon pisum) in these plants has no effect on ﬁne forceps in saline solution. The dissected organs were Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 1555 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Boulain et al. GBE soaked in RNA later (QIAGEN) immediately after dissection to Saliva and Salivary Gland Proteomics avoid RNA degradation. The dissected tissues were pooled in Saliva and Salivary Gland Collection several batches, and RNA was extracted by NucleoSpin RNA LSR1 aphids of mixed ages were reared on V. faba and ap- XS (Macherey-Nagel) and quantiﬁed. On average, RNA sam- proximately 2,000 aphids were installed on 12 perspex rings ples from 200 pairs of salivary glands or 20 alimentary tracts (radius 4.5 cm, height 5 cm), each containing 5 ml of a chem- that were dissected on the same day were pooled for one ically deﬁned diet formulation AP3 (Febvay et al. 1988)held replicate of RNA-seq experiment. Four replicates were pre- TM between two stretched sheets of Paraﬁlm . The aphids were pared by 4 days of dissection with 5–6 persons. reared at 18 C with 16 h day/8 h night photoperiod and the rRNA depletion, single stranded-RNA library preparation, diets were collected and replaced every 24 h. New V. faba- multiplexing, and sequencing were performed by Genewiz reared aphids were added to the diets to maintain aphid (New Jersey, USA). Sequencing was performed on the numbers. The daily collected diets were pooled and stored Illumina HiSeq2500 platform, in a 2 125 bp paired-end (PE) at 80 C for later use. Three independent replicates were conﬁguration in High Output mode (V4 chemistry). Each sam- produced by pooling the collected diet from three daily col- ple was sequenced on four different ﬂowcell lanes to avoid lections (approximately 150 ml). Pooled diets were concen- lane effect. In total, 471,933,074 reads were obtained for trated at 4 C in a Vivacell 250 Pressure Concentrator using eight samples. Raw data is available in NCBI Sequence a 5,000 molecular weight cut-off (MWCO) polyethersulfone Read Archive (https://trace.ncbi.nlm.nih.gov/Traces/sra/; last (PES) membrane. The volume of concentrate was further re- accessed May 22, 2018) with reference number SRP14110. duced by centrifugation at 3,400 g ina Vivaspin 6with a 5000 MWCO. Proteins from this ﬁnal concentrate (300 ml) Mapping and Differential Expression Analysis were puriﬁed using a 2D Clean-up Kit (GE HealthCare) fol- lowing the manufacturer’s instructions and the resulting pro- Gene expression of A. pisum salivary glands and alimentary tein pellet was suspended in 6 M urea, 2 M thiourea, 0.1 M tracts was analyzed using the Acyr_2.0 (GCF_000142985.2) TM Tris–HCl, pH 8.0, and quantiﬁed using the Qubit protein reference genome assembly and the NCBI Acyrthosiphon quantiﬁcation system (Invitrogen). Ten micrograms were pisum Annotation Release 102, both available at ftp://ftp. removed from each sample for protein digestion. ncbi.nlm.nih.gov/genomes, last accessed May 22, 2018. Adult LSR1 aphids (14–16-days old, reared on V. faba) The paired-end libraries were mapped on the reference were dissected in ice-cold saline and dissected SGs were im- genome using STAR v2.5.2 (Dobin et al. 2013) with mediately transferred to 60 ml PBS supplemented with Roche the following parameters: outFilterMultimapNmax¼ 5, TM cOmplete protease inhibitor cocktail (PIC; ﬁnal EDTA con- outFilterMismatchNmax¼ 3, alignIntronMin¼ 10, centration: 0.2 mM). Thirty pairs of SGs were pooled per rep- alignIntronMax¼ 50,000, alignMatesGapMax¼ 50,000. licate and homogenized with a disposable pestle. Sixty Fragment counts per genes were estimated by Subread microliter of 12 M urea, 4 M thiourea, and PIC was added featureCounts (Liao et al. 2014) using default parameters. and samples were homogenized further, centrifuged at Differential expression analysis between SGs and AT was 9,000 g for 5 min to pellet cellular debris and the superna- then conducted following the workﬂow proposed by Law tant was removed and quantiﬁed. One hundred microgram et al. (2016). The raw fragment counts were converted to of protein was removed and puriﬁed using the 2D Clean-up counts per million (CPM) using the edgeR (Robinson et al. Kit. The solubilized protein lysates were resuspended in 6 M 2010) R-implemented package (R-Core Team 2017). urea, 2 M thiourea, 0.1 M Tris–HCl, pH 8.0, and requantiﬁed. Expressed genes were ﬁltered based on a CPM>1in at Twenty micrograms were removed from each sample for pro- least three libraries among the eight analyzed libraries and tein digestion. CPMs were normalized by the edgeR TMM method for Normalization Factor calculation (Robinson and Oshlack Sample Preparation for Mass Spectrometry 2010). The mean–variance relationship of the log-CPM was estimated by the voom function (Law et al. 2014) For mass spectrometry, three independent biological repli- and incorporated in the empirical Bayes analysis from the cates were analyzed. Fifty-mM ammonium bicarbonate was limma package (Ritchie et al. 2015) to ﬁt linear models and added to each sample, and proteins were reduced with 0.5 M compare SG vs. AT tissues. Validation of the described dif- dithiothreitol (DTT) at 56 C for 20 min. The proteins were ferential expression analysis is presented in supplementary then alkylated with 0.55 M iodoacetamide (IAA) at room tem- ﬁg. S1, Supplementary Material online. In addition, normal- perature for 15 min, in the dark. One microliter of a 1% w/v ized fragments per kilobase million (FPKM) were calculated solution of Protease Max Surfactant Trypsin Enhancer with edgeR on the four salivary glands samples and the four (Promega) and 0.5 mg of Sequence Grade Trypsin (Promega) alimentary tracts samples separately to obtain whole tran- was added to obtain a protein: trypsin ratio of 40:1 and 80:1 scriptomes of each organ. Genes with FPKM> 0.5 in at for saliva and salivary glands, respectively. The protein/trypsin least three libraries per tissue were considered as expressed. mixture was incubated at 37 C for 18 h. Digestion was 1556 Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Evolution and Gene Family Expansions of Aphid Salivary Effectors GBE terminated by adding 1 ml of 100% triﬂuoroacetic acid (Sigma The data matrix was ﬁrst ﬁltered for the removal of contam- Aldrich) and incubation at room temperature for 5 min. inants and peptides identiﬁed by site. Label Free Samples were centrifuged for 10 min at 13,000 g and pu- Quantitation (LFQ) intensity values were log2 transformed riﬁed for mass spectrometry using either the ZipTip pipette and proteins not found in all three replicates in at least one procedure (Millipore) for saliva and C18 Spin Columns group were omitted from the analysis. A data-imputation (Pierce) for salivary glands, following the manufacturer’s step was conducted to replace missing values with values instructions. The eluted peptides were dried using a that simulate signals of low abundant proteins chosen ran- SpeedyVac concentrator (Thermo Scientiﬁc Savant domly from a distribution speciﬁed by a downshift of 1.8 DNA120) and resuspended in 2% v/v acetonitrile and times the mean standard deviation (SD) of all measured 0.05% v/v triﬂuoroacetic acid (TFA). Samples were sonicated values and a width of 0.3 times this SD. for 5 min to aid peptide resuspension followed by centrifu- gation for 5 min at 13,000 g. The supernatant was re- Gene Ontology and Secretion Prediction moved and used for mass spectrometry. The GO annotation was performed on the whole A. pisum proteome available on NCBI using Blast2GO v2.5.0 (Go¨tz et al. 2008). Associations were realized using a blastp Mass Spectrometry [BLASTþ v2.5.0 (Camacho et al. 2009)] search against the One microgram of each digested sample was loaded onto a nonredundant protein database (release 2017-2-4) with the QExactive (ThermoFisher Scientiﬁc) high-resolution accurate following parameters: e-value¼ 1e8, max_target_ mass spectrometer connected to a Dionex Ultimate 3000 seqs¼ 20, soft_making¼ false, and Interproscan v5.13.52.0 (RSLCnano) chromatography system. The peptides were sep- (Jones et al. 2014). Assigned GO terms for genes of interest arated by a 2–40% gradient of acetonitrile on a Biobasic C18 groups were categorized by molecular function (MF), biolog- PicofritTM column (100 mm length, 75 mm ID), using a 120- ical process (BP) and cellular component (CC) and GO enrich- min reverse-phase gradient at a ﬂow rate of 250 nl min ment was investigated using hypergeometric tests in R. The with a runtime of 50 and 130 min for saliva and salivary P-values were adjusted for multiple comparisons using glands, respectively. All data were acquired with the mass Bonferroni correction. The whole A. pisum proteome was spectrometer operating in automatic data dependent switch- analyzed with SignalP v3.0 and v4.1 (Dyrløv Bendtsen et al. ing mode. A full MS scan at 70,000 resolution and a range of 2004; Petersen et al. 2011) and SecretomeP v2.0 (Bendtsen 400–1,600 m/z was followed by an MS/MS scan at resolution et al. 2004) to characterize the presence of signal peptide or 17,500 and a range of 200–2,000 m/z, selecting the 15 most nonclassical secretion signal, respectively. Then, membrane intense ions prior to MS/MS. inserted domains were predicted using TMHMM v2.0 Protein identiﬁcation of MS/MS data was performed using (Krogh et al. 2001) for transmembrane domains and MaxQuant v1.5.6.5 (www.maxquant.org; last accessed May PredGPI (Pierleoni et al. 2008) for GPI anchors. Finally, the 22, 2018) following the general procedures and settings out- results were combined to deﬁne the list of A. pisum secreted lined in Hubner et al. (2010). The Andromeda search algo- proteins that have secretion signals and no membrane inser- rithm (Cox et al. 2011) implemented in the MaxQuant tion domains. software was used to correlate MS/MS data against the protein reference sequence set of A. pisum obtained from Orthology Analysis the NCBI (27,984 entries, May 2016) and a contaminant sequence set provided by MaxQuant. The following To determine groups of orthologs among the 17 genomes search parameters were used: ﬁrst search peptide toler- (supplementary table S1, Supplementary Material online), ance of 20 ppm, second search peptide tolerance 4.5 ppm we ﬁrst kept the longest protein isoforms from each species with cysteine carbamidomethylation as a ﬁxed modiﬁca- with an home-made perl script and then ran the tion and N-acetylation of protein and oxidation of methi- OrthoDB_soft_1.6 (Kriventseva et al. 2015) using standard onine as variable modiﬁcations and a maximum of two parameters. To establish the species phylogeny, 478 groups missed cleavage sites allowed. False Discovery Rates of conserved single-copy orthologs were extracted and their (FDR) were set to 1% for both peptides and proteins protein sequences aligned with MAFFT (Katoh et al. 2005). and the FDR was estimated following searches against a Alignments were concatenated and used to generate a target-decoy database. Peptides with minimum length of maximum likelihood tree with RAxML (Stamatakis 2006) seven amino acids were considered for identiﬁcation. with default parameters and 1,000 replicates, considering Tetranychus urticae as an outgroup. The phylogeny was then used to establish the level of orthology of each A. Data Analysis pisum gene. The enrichments of orthologous categories Perseus v.1.5.5.3 (www.maxquant.org, last accessed May 22, among data set were analyzed using hypergeometric test 2018) was used for data analysis, processing and visualization. implemented in R. Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 1557 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Boulain et al. GBE Evolution Analyses of Salivary Genes multigene families than expected by chance. For each catalog, a signiﬁcant effect is demonstrated if the number of genes Single Copy Genes that belong to multigene families that are represented by two Ortholog groups that were represented by a single gene in the or more copies in the catalogue lies above the 95% conﬁ- A. pisum genome and in at least one of the other Aphididae dence interval (CI) of the expectation. 95% CI was computed genomes were extracted. If alternative transcripts of a gene by randomly sampling the number of genes contained in a were present, the longest coding sequence (CDS) for each catalog (740 and 3,603 genes) from the list of 18,603 genes species was used. Pairs of orthologous CDS (one CDS per and counting the number of gene-family members (that be- species: A. pisum/A. glycines, A. pisum/D. noxia, A. pisum/ long to a multigene family represented by two more copies in M. persicae) were aligned using MAFFT (Katoh et al. 2005) the sample) in this random sample. This step was repeated within the Hot algorithm from GUIDANCE2 (Sela et al. 2015) 10,000 times. to produce reliable codon-based alignments without poorly To test for positive selection on certain sites and branches aligned regions. From these codon alignments, pairwise dN/ (selection acting on foreground branch compared with back- dS values were calculated with PAML v4.8 using the YN00 ground branches) of the three selected gene families program (Yang and Nielsen 2000; Yang 2007). We removed (Cysteine-Rich Protein, Angiotensin-Converting-Enzyme-like ortholog pairs with dS> 2or dN> 0.5toavoid mutational and Aminopeptidase-N families), a cleaned alignment and a saturation, as well as pairs with dS¼ 0. dN/dS of different maximum-likelihood tree were generated as described above, gene categories within each species pairs were compared and were fed to the codeml branch-site (BS) model (Yang and using a Kruskal–Wallis test, and Nemenyi-Tests for multiple Nielsen 2002). The BS model classiﬁed the sites in four cate- comparisons were realized with the R package PMCMR gories: class 0 where sites were under negative selection (Pohlert 2014). (x0< 1) on both foreground and background branches, class Orthologs of characterized salivary effectors were searched 1 where sites evolved under neutral evolution (x1¼ 1) on in all Aphididae sequences available in Genbank database and both foreground and background branches, class 2a where in 454-sequenced contigs of Rhopalosiphum padi and the sites were under positive selection (x2 1) on the fore- Schizaphis graminum assembled for this study (see supple- ground branch and under negative selection (x0< 1) on mentary material and methods S1, Supplementary Material background branches, and class 2b where the sites evolved online), in order to compute more accurate evolutionary rates under positive selection (x2 1) on the foreground branch and to test selection models. For this task, reciprocal blast and under neutral evolution (x1¼ 1) on background searches were runwithblastn(BLASTþ v2.5.0) using A. branches. The null model in which the foreground branch pisum, M. persicae, A. glycines and D. noxia CDS sequences may have different proportions of sites under neutral evolu- against the different databases (e-value< 10 ). Then, for tion (dN/dS of 1), and the alternative model, in which the each obtained ortholog group, a codon-based alignments foreground branch may have sites under positive selection, was generated with PRANK (Lo ¨ ytynoja and Goldman 2005) were applied to all branches of each gene family. The likeli- and cleaned from poorly aligned regions using GUIDANCE2, hood of the selection model was computed as described in specifying a minimum alignment quality threshold of 0.93. previous section and P-values were corrected for multiple test- The cleaned alignment was used to generate a maximum ing comparisons as described by Anisimova and Yang (2007). likelihood phylogenetic tree with RAxML and both alignment In case the selection model was retained, the posterior prob- andtree servedtoestimate dN/dS with codeml [implemented ability of particular sites evolving under positive selection was in PAML v4.8 (Yang 2007)] under different models (Yang et al. computed. 2000; Yang and Nielsen 2000). The null models (M0¼ one All phylogenetic trees shown in ﬁgures were designed in ratio, M1¼ neutral, M7 ¼ b) were compared with alternative iTOL (Letunic and Bork 2007). models (M2¼ selection, M8 ¼ bþ x) using the likelihood ratio test (LRT), which compares twice the difference in log Results likelihood to a v distribution. Finally, if selection models were found more likely, the Bayes Empirical Bayes (BEB) procedure Generation of the A. pisum Salivary Effector Candidate was used to compute the posterior probability of evolution Gene Sets under positive selection for each site (Yang et al. 2005). We conducted transcriptome analyses of A. pisum salivary glands (SGs) to generate two catalogues of salivary effector Gene Families candidates: one that considers up-regulation of their expres- Families containing genes expressed in A. pisum salivary sion in SGs and another that does not. This approach assumes glands and their orthologs in arthropod species were that the vast majority of salivary proteins secreted into plants extracted from OrthoDB groups. Then, we examined if the are expressedinaphidSGs (Mutti et al. 2008; Carolan et al. two salivary gene catalogs contain more members of 2011; Naessens et al. 2015) and that most salivary effectors 1558 Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Evolution and Gene Family Expansions of Aphid Salivary Effectors GBE AB FIG.1.—Pipelines used to establish sets of candidate Acyrthosiphon pisum salivary effector genes expressed and up-regulated in salivary glands (A). Pipeline used to identify proteins from saliva injected in artiﬁcial diet and composition of the secreted proteins (B). are highly expressed in SGs compared with other organs compared between SGs and AT. After ﬁltering and normali- (Wang et al. 2015a, 2015b). Up-regulation of expression in zation steps, 12,378 protein coding genes from SGs and AT SGs was determined in comparison to expression levels in the were retained. Among these genes, 1,989 genes were up- alimentary tract (AT), which we chose due to ease of isolation regulated in SG tissues, of which 740 genes were predicted to and RNA extraction. RNA-seq by Illumina technology was encode secreted proteins. These constituted the candidate conducted on dissected SG and AT tissues of A. pisum.In SG-up effector set (ﬁg. 1A and supplementary table S2, SGs, 12,040 genes passed the cut-off value for gene expres- Supplementary Material online). All were found in the sion (ﬁg. 1A) and encode proteins that may or may not be 3,603 SG-expressed effector set except one (gene secreted in saliva. A protein secretion prediction pipeline was LOC100574271), which was expressed at too low level to applied to the global gene set of A. pisum (N¼ 18,601) using pass the criterion used to deﬁne the SG-expressed effector SignalP v3 or v4.1 (Dyrløv Bendtsen et al. 2004; Petersen et al. set (FPKM 0.5) although it passed the criterion used to de- 2011), as well as SecretomeP (Bendtsen et al. 2004), which ﬁne the SG-up effector set (CPM> 1). predicts noncanonical secretion signals. Genes encoding transmembrane domains [predicted by TMHMM 2.0 (Krogh Identiﬁcation of Aphid Proteins in Artiﬁcial Diets and et al. 2001)] or GPI-anchors [predicted with PredGPI (Pierleoni Salivary Glands et al. 2008)] were considered as not secreted, leaving 3,603 encoded proteins predicted to be secreted. These constituted Proteins from artiﬁcial diets fed upon by aphids were analyzed the candidate SG-expressed effector set (ﬁg. 1A and supple- by proteomics-based mass spectrometry (MS). Fifty-one pro- mentary table S2, Supplementary Material online). teins were supported by more than one peptide detected by To create the candidate SG-up-regulated (or SG-up) effec- MS in at least two out of the three replicates or by one peptide tor set, the expression levels of individual RNAs were in all three replicates of saliva samples from artiﬁcial diets Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 1559 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Boulain et al. GBE (ﬁg. 1B and supplementary table S3, Supplementary Material estimating the ratio of nonsynonymous to synonymous online). Of these, 35 proteins (68.62%) were included in the substitutions (dN/dS or x) in three categories: SG-up ef- SG-up effector set, two belonged to the SG-expressed effec- fector set, SG-expressed effector set and the remaining tor set, and 14 were not predicted as secreted proteins and genes (ﬁg. 3). Only single-copy gene pairs of orthologs were encoded by genes up-regulated in SGs (11 genes) or between A. pisum and the three other Aphididae species expressed in SGs (three genes) (ﬁg. 1B and supplementary were included in the orthology search, and pairwise evo- table S2, Supplementary Material online). Five out of the 11 lutionary rates were analyzed (A. pisum/A. glycines, A. SG-up-regulated genes seemed to be truncated in the ge- pisum/D. noxia, A. pisum/M. persicae). On average, genes nome sequence and were predicted to be not secreted. of SG-up effector set showed higher dN/dS ratios than the In all, 1,837 proteins (supported by more than one peptide genes of other sets (ﬁg. 3). The same analysis comparing in at least two of the three replicates or by one peptide in all AT-up, AT-expressed and remaining genes revealed that the three replicates) encoded by 1,809 unique genes were the observed overall higher dN/dS ratio was speciﬁc to the detected directly from dissected SGs (supplementary table genes of SG-up effector set (supplementary ﬁg. S3B, S3, Supplementary Material online), among which 205 pro- Supplementary Material online). teins belonged to the SG-up effector set and 447 proteins to The 11 single-copy salivary effector genes that were previ- the SG-expressed effector set. Signal intensity of proteins in ously functionally characterized and shown to be involved in the SGs signiﬁcantly correlated (Pearson’s r¼ 0.5446, plant-aphid interactions were examined (table 1). All these P< 0.001) with the transcription level of corresponding genes single copy salivary effectors or their A. pisum orthologs in the organ although such correlation was not observed for were found in the SG-up effector set with the exception of the proteins detected in the artiﬁcial diets (supplementary ﬁg. Mif1 (Naessens et al. 2015), which was not up-regulated in S2, Supplementary Material online). SGs and not predicted to be secreted in our study (table 1). Most of these genes were highly expressed, and genes C002, Mp1, Ap25, Me23, Me10, Mp55, Shp,and Mp2 were Aphididae Speciﬁc Genes Are Enriched in the SG-Up- Aphididae- or Aphidomorpha-speciﬁc. Regulated Effector Set To examine the evolutionary rates of these 11 effectors, we To assess the phylogenetic distribution of genes constituting retrieved the orthologs from SG transcripts of two additional the SG-up and SG-expressed effector sets, an analysis of aphid species, Schizaphis graminum and Rhopalosiphum padi orthology was conducted between A. pisum and16other (supplementary material and methods S1 and table S4, arthropods whose genome sequences were available, and Supplementary Material online). In addition, orthologs from which cover a wide range of divergence from A. pisum other aphids were retrieved from public databases (supple- (ﬁg. 2A). Relative to the whole set of genes annotated in mentary table S4, Supplementary Material online) and used the A. pisum genome, the SG-expressed effector set was to estimate dN/dS ratios under different site-models of selec- signiﬁcantly enriched in highly conserved genes found tion. The global dN/dS for each gene showed a range from across the arthropod or insect clades (hypergeometric 0.06298 (Armet) to 0.64653 (Mp1)(table 1 and supplementary test, P< 0.01). The set was also slightly enriched in genes table S5, Supplementary Material online). Differences were found only in Aphidomorpha, whereas the proportion of A. found between insect- or arthropod-conserved effector genes pisum-speciﬁc genes was signiﬁcantly reduced (hypergeo- like Armet, Mif1 and Mp10, which seem to have evolved under metric test, P< 0.001). In contrast, in the SG-up effector strong negative selection (dN/dS< 0.12), and the aphid- set, the proportion of genes found only in Aphididae was speciﬁc, fast-evolving effector genes (dN/dS> 0.40, in the signiﬁcantly higher than in the whole gene set, and highly top 5% of values obtained from pairwise comparisons). conserved genes found across arthropods were less fre- Despite these faster evolutionary rates, the examination of quent (hypergeometric test, P< 0.001) (ﬁg. 2B). In our com- null and alternative models of codon substitutions [M1 vs. parison, this pattern was speciﬁc to the SG-up effector set M2 and M7 vs. M8] revealed signatures of positive selection as the genes expressed or up-regulated in ATs showed sig- only for genes Mp1 and Me10 (supplementary table S5, niﬁcant enrichment of highly conserved genes found across Supplementary Material online). In Mp1 and Me10,sites arthropod and insect clades (supplementary ﬁg. S3A, detected as under positive selection were mainly found in the Supplementary Material online). region coding for the mature protein injected into the host plant (supplementary ﬁg. S4, Supplementary Material online). Evolutionary rates vary among of A. pisum Salivary Candidate Effectors Episodic positive selection occurred in salivary expanded The enrichment of Aphididae-speciﬁc genes in the SG-up gene families in the A. pisum Lineage effector set may be associated with rapid molecular evo- lution related to aphid-speciﬁc gene functions. This hy- Evolutionary novelties are usually brought by gene duplication pothesis was tested for the SG-up effector set by followed by diversiﬁcation in functions (Conant and Wolfe 1560 Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Evolution and Gene Family Expansions of Aphid Salivary Effectors GBE Tetranychus urticae Nasonia vitripennis Apis mellifera Tribolium castaneum Arthropoda Bombyx mori Manduca sexta Anopheles gambiae Drosophila melanogaster Insecta Mayetiola destructor Nilaparvata lugens Rhodnius prolixus Cimex lectularius Hemiptera Daktulosphaira vitifoliae 100 Aphis glycines Aphidomorpha 0.1 Acyrthosiphon pisum Aphididae Myzus persicae Diuraphis noxia SG-up-regulated set (740) *** *** Arthropoda Insecta Hemiptera SG-expressed set Aphidomorpha (3603) *** * * *** Aphididae A. pisum specific Whole gene set (18,601) 0 20 40 60 80 100 Percentage of genes FIG.2.—Orthology proﬁle of Acyrthosiphon pisum SG effector candidate gene sets (478 single-copy ortholog groups). (A) Phylogeny of the 17 arthropod species analyzed to determine ortholog groups. The colored squares indicate the levels of orthology, as indicated on the bottom right-hand legend. (B) Proportions of the different orthology levels among genes of the salivary effector sets and the A. pisum genome. Stars indicate the signiﬁcance of differences in the proportion of genes presenting a given orthology level between a given effector set and the whole gene set (hypergeometric test): *P< 0.05, ***P< 0.001. 2008; Kondrashov 2012). Interestingly, duplicated genes are in proteolysis, two functions for which the SG-up effector more frequently observed in SG-up and SG-expressed effector set was enriched (supplementary ﬁg. S5, Supplementary sets than expected. Indeed, the observed number of genes Material online). that are represented by two or more members of gene family in the two catalogs always lies above the 95% conﬁdence The Cysteine-Rich Protein (CRP) Family interval (CI) (SG-up: 80 genes, 95% CI ¼ [23, 55] and SG- expressed: 654 genes, 95% CI ¼ [442, 540]). Three multi- The CRP Aphididae speciﬁc gene family was originally de- gene families in particular, a cystein-rich protein family scribed as a family of 12 genes in A. pisum (Guo et al. (CRP), the Angiotensin-Converting-Enzyme-like family 2014). Here, 15 gene members were identiﬁed from genomic (ACE) and the Aminopeptidase-N (apN) family, represented data, and encode proteins of <200amino acidswith14 46.9% of the proteins detected in A. pisum saliva (supple- highly conserved cysteine residues. The family is expanded mentary table S6, Supplementary Material online). Gene in A. pisum with 15 copies compared with the other aphid ontology analysis showed that the ACE and apN gene fam- species D. noxia, A. glycines and M. persicae, which have 6, ilies belong to the metallopeptidase family and are involved 4, and 1 copies, respectively. Among the 15 A. pisum gene Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 1561 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Boulain et al. GBE copies, three were not detected in either SG or AT tissues Others SG-expressed candidates by RNAseq or protein mass spectrometry. Eleven genes SG-up-regulated candidates were up-regulated in SGs, nine encode proteins that were predicted to be secreted (supplementary table S6, a b Supplementary Material online), two of which (CRP-12 a and 14) were found in saliva. The branch-site (BS) model of codon substitutions (Yang and Nielsen 2002) detected positive selection in ﬁve branches of the Aphididae CRP family tree (ﬁg. 4A and B): the branch leading to D. noxia gene CRP-5 (branch #1), the ancestral branches #2, #4, and those leading to A. pisum CRP-2 (branch #3) and CRP-12 (branch #5). The Bayes Empirical Bayes (BEB) pos- A. pisum/A. glycines A. pisum/D. noxia A. pisum/M. persicae terior probability (Yang et al. 2005) revealed sites under positive selection (BEB above 0.75) in each aforemen- FIG.3.—Pairwise evolutionary rates (dN/dS) of single-copy Aphididae tioned branch (supplementary ﬁg. S6A, Supplementary genes orthologous to Acyrthosiphon pisum effector candidates. Here, SG- Material online), which are scattered across the protein expressed effector set does not contain SG-up effector set. For the three (supplementary ﬁg. S6B, Supplementary Material online). species pairs indicated on the X axis 5297, 5379, and 5562 single copy The functional importance of these sites is unknown, ortholog pairs were retained, respectively. Letters above boxes denote signiﬁcant differences determined by multiple Kruskal–Wallis test within due to the absence of known domains within the proteins. each species pair (A. pisum/Aphis glycines: H¼ 33.944, 2 d.f., P< 0.001; Although positive selection seems to occur in some A. A. pisum/Diuraphis noxia: H¼ 51.17, 2 d.f., P< 0.001;A.pisum/Myzus pisum CRP copies, no site under positive selection persicae: H¼ 39.373, 2 d.f., P< 0.001). was detected in CRP-13, the most highly expressed A. pisum CRP. Table 1 Aphid Single-Copy Salivary Effector Genes Characterized for their Role in the Interaction with Host Plants Gene Acyrthosiphon Effector set Expression Orthology x dN/dS Aphid Phenotype References* pisum Gene Rank Level Mp1 LOC100165393 SG-up 19/740 Aphididae 0.64653 Expression promotes fecundity [3, 4, 6, 10, 15] ACYPI006346 C002 LOC100167863 SG-up 12/740 Aphididae 0.62307 Expression promotes fecundity [1, 2, 3, 4, 6, 13] ACYPI008617 Silencing reduces survival/fecundity Ap25 LOC100169287 SG-up 57/740 Aphididae 0.56219 Expression promotes fecundity [14] ACYPI009919 Me23 LOC100161198 SG-up 109/740 Aphididae 0.53025 Expression promotes fecundity [5] ACYPI002439 Me10 LOC100167427 SG-up 5/740 Aphididae 0.51787 Expression promotes fecundity [5] ACYPI008224 Shp LOC100169243 SG-up 1/740 Aphidomorpha 0.48863 Silencing reduces survival/fecundity [8, 12] ACYPI009881 Mp55 LOC100569515 SG-up 18/740 Aphididae 0.42972 Expression promotes fecundity [7] ACYPI33755 Silencing reduces survival/fecundity Mp2 LOC100160479 SG-up 107/740 Aphididae 0.40948 Expression promotes fecundity [6] ACYPI001774 Silencing reduces survival/fecundity Mp10 LOC100145855 SG-up 303/740 Insect 0.12212 Expression reduces fecundity [3] ACYPI000097 Mif1 LOC100161225 None 4700/12040 Arthropod 0.09454 Silencing reduces fecundity [9] ACYPI002465 Armet LOC100167188 SG-up 122/740 Arthropod 0.06298 Silencing reduces survival [11] ACYPI008001 Expression ranks show the rank of effector SG expression level (based on Log[FPKM]) within the SG-up effector set except for Mif1 which does not belong to any effector sets (global SG expression scale). x values represent the evolutionary rates within closely related Aphididae species. 1 2 3 4 5 6 7 8 * Mutti et al. 2006; Mutti et al. 2008; Bosetal. 2010; Pitino et al. 2011; Atamian et al. 2013; Pitino and Hogenhout 2013; Elzinga et al. 2014; Abdellatef et al. 2015; 9 10 11 12 13 14 15 Naessens et al. 2015; Pan et al. 2015; Wang et al. 2015a; Will and Vilcinskas 2015; Zhang et al. 2015; Guy et al. 2016; Rodriguez et al. 2017. 1562 Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 dN/dS 0.0 0.5 1.0 1.5 2.0 Evolution and Gene Family Expansions of Aphid Salivary Effectors GBE FIG.4.—Evolutionary analysis of the Aphididae CRP gene family. (A) Codon-based maximum-likelihood tree used to compute dN/dS under the branch- site (BS) model of codon substitution. Numbers below branches indicate bootstrap support (percentage of the 1,000 replicates) when >40. DNOX, MPER, AGLY, and APIS indicate sequences from Diuraphis noxia, Myzus persicae, Aphis glycines,and Acyrthosiphon pisum, respectively. Branches where positive selection affected certain sites are in bold. (B) Estimated evolutionary parameters for these branches (numbered as on the tree). x0, x1, and x2 indicate the average dN/dS ratiofor sites assignedtoclass 0(x < 1, negative selection), class 1 (x ¼ 1, neutral evolution) and class 2 (x > 1, positive selection), respectively, in the branch, and n.c. indicates insufﬁcient dS to compute x. p0, p2a, and p2b show proportions of sites in classes 0, 2a, and 2b, respectively, for each branch (Yang and Nielsen 2002). genes were up-regulated in A. pisum SG tissue, including The Angiotensin-Converting-Enzyme-like (ACE) Gene three members (ACE1, ACE2,and ACE5) of clade 4 encoding Family predictedsecretedACEs. ACE1, 5 and10were detectedin In contrast to the Aphididae-speciﬁc CRP gene family, the saliva, although the latter lacks an encoded signal peptide, ACE gene family is found among arthropods (ﬁg. 5 and sup- indicating errors in gene annotation or secretion prediction plementary table S6, Supplementary Material online). A (supplementary table S6, Supplementary Material online). It is search for the ortholog groups yielded ten members in the notable that the other aphid genomes comprise just one ACE A. pisum genome [three were previously identiﬁed by (Wang of clade 4 (ﬁg. 6A). All clade-4 ACE proteins except ACE10 et al. 2015b)], which are distributed in four clades along with are determined as functional peptidases, based on the pres- orthologs from other insect species (ﬁg. 5 and supplementary ence of the M2 peptidase HEXXH domain (IPR001548) table S6, Supplementary Material online). Seven of these (ﬁg. 6A). Positive selection was inferred in three branches of Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 1563 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 79 AMEL XP_006560388.1 RPRO RPRC009977-PA NVIT NV12891-PA CLEC XP_014260341.1 TCAS XP_974768.1 NLUG NLU015818.1 CLEC XP_014240472.1 APIS XP_008188413.1 ACE8 RPRO RPRC003668-PA APIS XP_001947791.1 ACE9 NLUG NLU018682.1 MPER G006 v1.0 000174470 APIS XP_001949396.4 ACE3 DNOX XP_015365527.1 DNOX XP_015374702.1 APIS XP_001943158.2 ACE2 MPER G006 v1.0 000139710 AGLY AG001409-PA APIS XP_016663856.1 ACE7 AGLY AG001407-PA AGLY AG006483-PA DVIT DV3001051-PA TURT tetur25g01290.1 DVIT DV3001049-PA NLUG NLU015034.2 RPRO RPRC007881-PA AGAM AGAP009756-PA CLEC XP_014258853.1 AGAM AGAP009755-PA RPRO RPRC001697-PA AGAM AGAP009754-PA Boulain et al. GBE SG up-regulated Undifferentiated SG-expressed effector set SG-up effector set ACE Detected in saliva clade 1 ACE clade 2 ACE clade 4 ACE clade 3 0.5 FIG.5.—Protein-based maximum likelihood phylogeny of the ACE gene family among insect species, with bootstrap support indicated next to branches (percentage of the 1,000 replicates) when >40. Species abbreviations: AGAM, Anopheles gambiae;AGLY, Aphis glycines;AMEL, Apis mellifera; APIS, Acyrthosiphon pisum;BMOR, Bombyx mori;CLEC, Cimex lectularius;DMEL, Drosophila melanogaster; DNOX, Diuraphis noxia;DVIT, Daktulosphaira vitifoliae;MDES, Mayetiola destructor;MPER, Myzus persicae; MSEX, Manduca sexta;NLUG, Nilaparvata lugens;NVIT, Nasonia vitripennis;RPRO, Rhodnius prolixus;TCAS, Tribolium castaneum;TURT, Tetranychus urticae. this clade (ﬁg. 6B and C): the ancestral branch (#1) of A. pisum The Aminopeptidase-N (apN) Gene Family ACE 1 and ACE5,the A. pisum ACE1 gene (#2), and the ancestral branch #3. Sites detected to evolve under positive The aminopeptidase-N (apN) gene family was the most rep- selection are scattered across the protein (supplementary ﬁg. resented in A. pisum saliva proteome with 18 apN proteins S7A and B, Supplementary Material online). detected, accounting for 36.7% of the 51 saliva proteins 1564 Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 NLUG NLU015031.1 AGAM AGAP009752-PA BMOR XP_012550154.1 AGAM AGAP009751-PA MSEX Msex2.07574-RB DMEL FBpp0080129 MSEX Msex2.07570-RB DMEL FBpp0079297 AGAM AGAP009757-PA MDES Mdes006926-RB DMEL FBpp0099733 BMOR Gene005248 TCAS NP_001164243.1 MSEX Msex2.07568-RG NVIT NV13720-PA MSEX Msex2.07569-RA AMEL XP_392082.3 BMOR Gene005250 DNOX XP_015365902.1 TCAS XP_008198828.2 APIS XP_001948448.3 ACE6 DNOX XP_015364631.1 MPER G006 v1.0 000080860 MPER G006 v1.0 000019250 AGLY AG006015-PA AGLY AG009522-PA DVIT DV3016127-PA APIS NP_001129384.1 ACE4 MDES Mdes005391-RA APIS XP_003242883.2 ACE1 AGAM AGAP007982-PA APIS XP_008188393.1 ACE5 MSEX Msex2.02138-RB APIS XP_001951640.3 ACE10 BMOR Gene000780 DVIT DV3004498-PA Evolution and Gene Family Expansions of Aphid Salivary Effectors GBE FIG.6.—Evolutionary analysis of clade 4 of the Aphididae ACE gene family. (A) Protein structure of the family members with protease active sites shown as blue lines. (B) Codon-based maximum likelihood tree used to compute dN/dS under the branch-site (BS) model of codon substitution. Numbers below branches indicate bootstrap support (percentage of the 1,000 replicates) when >40. DNOX, MPER, AGLY, and APIS indicate sequences from Diuraphis noxia, Myzus persicae, Aphis glycines and Acyrthosiphon pisum, respectively. Branches where positive selection was detected at speciﬁc sites are shown in bold. (C) Estimated evolutionary parameters for these branches (numbered as on the tree). See ﬁg. 4B for the deﬁnition of terms used. identiﬁed. Another apN protein not detected here was iden- followed by a M1 peptidase domain (IPR014782) and an tiﬁed in saliva by Carolan et al. (2009) (supplementary table ERAP1-C domain (IPR024571). However, some have lost S6, Supplementary Material online). Forty-seven apN genes one of the two domains, and the vast majority (apart from were identiﬁedinthe A. pisum genome, 27 of them belonged D. noxia apN-3, A. pisum apN-1 and A. pisum apN-5) pos- to ortholog groups present among insects, and the remaining sesses the HEXXH active site ensuring the peptidase function 20 did not show orthologs in other species (supplementary (ﬁg. 8A). table S6, Supplementary Material online). All 20 A. pisum- Traces of positive selection were detected by the branch- speciﬁc apN genes were up-regulated in SG tissue, ten en- site model in seven branches of apN clade 4 (ﬁg. 8B and C): code proteins with predicted secretion, four of which were the branch leading to D. noxia gene apN-1 (branch #1), the detected in saliva. Four of the ten proteins not predicted as ancestral branches #3, #5, #7, and those leading to A. pisum secreted were also detected in saliva, indicating incorrect apN-3 (#2), apN-8 (#4), and apN-12 (#6) genes. Along these gene annotation and/or secretion prediction. The 27 apN branches, sites under positive selection (with BEB> 0.75) genes with orthologs in other species were clustered in eight were mainly located in the M1 peptidase and ERAP1-C distinct clades (ﬁg. 7). The A. pisum genes from clades 1, 2, 3, domains, and in the uncharacterized part between the signal 4, and 7 were up-regulated in SGs with the exception of only peptide sequence and M1 peptidase domain for branches #1, three members. Genes from clades 5, 6, and 8 were not up- #2, and #3 (supplementary ﬁgs. S8 and S9, Supplementary regulated in SGs. Apart from the one A. pisum gene that Material online). belongs to clade 1, all the other SG-up candidate effector genes were found in clade 4 (ﬁg. 7). Clade 4 is remarkable for the presence of 14 closely related Discussion gene copies in A. pisum and only a few copies from other New Salivary Candidate Effector Gene Sets insect species, including aphids (ﬁg. 7). Eleven A. pisum copies As evidence accumulates that aphid salivary effectors play were SG-up-regulated and detected in saliva, but two encode important roles in plant-aphid interactions (Elzinga and proteins that were not predicted as secreted (ﬁg. 7). Genes of Jander 2013), a comprehensive identiﬁcation of aphid salivary this clade generally encode proteins with a signal peptide Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 1565 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 AGLY AG010842-PA APIS NP_001119606.1 DVIT DV3006347-PA APIS XP_001944286.2 MPER G006 v1.0 000108970 MPER G006 v1.0 000075720 AGLY AG000258-PA DNOX XP_015369492.1 DNOX XP_015371746.1 AGLY AG006939-PA APIS XP_008179876.1 DNOX XP_015378862.1 DVIT DV3004461-PA MPER G006 v1.0 000125280 CLEC XP_014254770.1 AGLY AG006771-PA RPRO RPRC002589-PA DNOX XP_015378863.1 NLUG NLU013342.1 MPER G006 v1.0 000125260 AMEL XP_392602.4 AGLY AG006770-PA NVIT NV10578-PA MPER G006 v1.0 000125270 BMOR Gene004818 DNOX XP_015378898.1 MSEX Msex2.05162-RA APIS XP_008179772.1 MDES Mdes013476-RA MPER G006 v1.0 000125300 AGAM AGAP000129-PA DNOX XP_015377443.1 DMEL FBpp0110455 AG006773-PA TCAS XP_008190425.1 DVIT DV3003804-PA MPER G006 v1.0 000118390 APIS XP_001946754.2 DNOX XP_015364307.1 APIS XP_001948442.2 APIS XP_016663568.1 MPER G006 v1.0 000144060 MPER G006 v1.0 000045750 DNOX XP_015378864.1 DNOX XP_015379576.1 MPER G006 v1.0 000144070 AGLY AG002007-PA APIS XP_001948350.2 100 86 DVIT DV3000262-PA DNOX XP_015378865 1 NLUG NLU008387.1 AGLY AG006769-PA NLUG NLU012822.1 NLUG NLU025909.2 NLUG NLU028309.1 APIS XP_001942841.2 CLEC XP_014243033.1 APIS XP_016662011.1 RPRO RPRC003230-PA APIS XP_016658913.1 BMOR Gene006567 APIS XP_016656075.1 MSEX Msex2.02406-RA APIS XP_003240500.1 TCAS XP_008200193.1 APIS XP_008189748.1 MDES Mdes000745-RA APIS XP_008179628.1 AGAM AGAP005728-PA APIS XP_016656878.1 DMEL FBpp0072583 APIS XP_008189321.2 AMEL XP_394245.2 Boulain et al. GBE SG up-regulated apN clade 2 AT up-regulated apN Undifferentiated apN clade 1 clade 3 SG-expressed effector set SG-up effector set Detected in saliva apN clade 4 apN clade 8 apN clade 7 apN clade 5 apN 0.5 clade 6 FIG.7.—Protein-based maximum likelihood phylogeny of the apN gene family among insect species, with bootstrap support indicated next to branches (percentage of the 1,000 replicates) when >40. Species abbreviations are the same as on ﬁg. 5. genes and analysis of their evolutionary patterns are essential sequencing. The salivary gland secretome established for A. to shed light on the process of aphid adaptation to speciﬁc pisum by Carolan et al. (2011) contains 324 genes, of which host plants. Here, we have assembled new catalogues of can- 95 are found in our SG-up effector set and 93 others are didate A. pisum salivary effector genes based primarily on contained in the SG-expressed effector set. The remaining comparative transcriptional analysis of SG and AT genes. 120 encompass 111 SG-expressed genes encoding proteins We identiﬁed 3,603 genes expressed in SGs and potentially considered as not secreted in our analysis, and 9 genes that secreted in saliva. A subset of 740 genes have their expression were absent from our source transcriptome. Finally, 16 up-regulated in SGs compared with AT. These sets combined genes did not exist in the updated gene annotation of the are substantially larger than the previously established candi- A. pisum genome used here (NCBI Acyrthosiphon pisum date effector catalogue (Carolan et al. 2011), reﬂecting the Annotation Release 102). The aphid line used here (LSR1) greater sensitivity of Illumina’s RNA-seq compared with EST was collected from Medicago sativa (The International 1566 Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 APIS XP_001948848.3 NVIT NV18310-PA APIS XP_008182457.2 NVIT NV16112-PA APIS XP_001950011.2 TURT tetur01g06430.1 DNOX XP_015367111.1 APIS XP_003242073.1 DNOX XP_015377346.1 AG010517-PA MPER G006 v1.0 000124450 MPER G006 v1.0 000003470 APIS XP_001950046.4 DNOX XP_015367797.1 APIS XP_001944764.2 DVIT DV3013700-PA DNOX XP_015375929.1 RPRO RPRC005596-PA DVIT DV3000933-PA CLEC XP_014248442.1 NLUG NLU028994.1 NLUG NLU011945.1 CLEC XP_014240291.1 TCAS XP_008200585.1 TCAS XP_008200837.1 AMEL XP_006571773.1 DMEL FBpp0084770 NVIT NV16959-PA DMEL FBpp0084767 TURT tetur02g06880.1 DMEL FBpp0084772 TURT tetur02g06960.1 MDES Mdes011987-RA BMOR Gene006183 AGAM AGAP000885-PA MPER G006 v1.0 000095660 BMOR Gene014807 MPER G006 v1.0 000095670 MPER G006 v1.0 000190520 APIS NP_001192196.1 APIS XP_016660197.1 MPER G006 v1.0 000095650 DNOX XP_015376659.1 DNOX XP_015379483.1 AGLY AG008289-PA APIS XP_003240091.1 DVIT DV3007457-PA MPER G006 v1.0 000095620 CLEC XP_014255314.1 AGLY AG006946-PA RPRO RPRC012350-PA DVIT DV3007809-PA NLUG NLU018374.3 CLEC XP_014246646.1 NVIT NV10584-PA CLEC XP_014246802.1 AMEL XP_006565537.1 NLUG NLU002176.1 AGAM AGAP013001-PA AMEL XP_006571893.1 DMEL FBpp0084751 NVIT NV14979-PA MDES Mdes001404-RA AGAM AGAP003077-PB MSEX Msex2.10159-RA DMEL FBpp0288476 TCAS XP_008200957.1 DMEL FBpp0082211 TURT tetur21g03280.1 DMEL FBpp0082209 MPER G006 v1.0 000178460 TCAS XP_008191642.1 APIS XP_008185292.1 MSEX Msex2.11506-RC DNOX XP_015363488.1 BMOR Gene007290 Evolution and Gene Family Expansions of Aphid Salivary Effectors GBE AB FIG.8.—Evolutionary analysis of the Aphididae clade 4 of the apN gene family. (A) Protein structure of the family members. (B) Codon-based maximum likelihood tree used to compute dN/dS under the branch-site (BS) model of codon substitution. Numbers below branches indicate bootstrap support (percentage of the 1,000 replicates) when>40. DNOX, MPER, AGLY, and APIS indicate sequences from Diuraphis noxia, Myzus persicae, Aphis glycines,and Acyrthosiphon pisum, respectively. The branches where positive selection was detected at certain sites are in bold. (C) Estimated parameters of these branches (numbered as on the tree). See ﬁg. 4B for the deﬁnition of terms used. Aphid Genomics Consortium 2010) and feeds well on the responses against aphids (Naessens et al. 2015). In addition, universal host plant of A. pisum, Vicia faba (our unpublished some effector proteins could be produced in tissues other data). We used V. faba to rear the aphids for RNA-seq ex- than SGs and injected into plants. Proteins can move from periment because of the ease of producing large quantity of the hemocoel of aphids to SGs to be secreted through saliva. synchronized aphids on this plant. It was also because pre- However, only the chaperonin GroEL, which is produced by vious study revealed that A. pisum reared on two suitable the aphid obligate endosymbiont in bacteriocytes, is known to legume plants show very similar expression patterns of sal- be injected into plant tissue through this pathway (Chaudhary ivary effector candidates (Eyres et al. 2016). There is a pos- et al. 2014). sibility that our catalogues missed the salivary genes that As 46 proteins out of the 51 (90%) detected in saliva were expressed exclusively when the aphid was feeding through proteomics were encoded by SG-up genes, most on M. sativa, but we assume that such genes are rare. of A. pisum salivary proteins seem encoded by genes whose Comparing expression levels between SGs and AT was expression is up-regulated in SGs. The lack of signiﬁcant fundamental in restricting the list of candidate effectors, correlation between the signal intensities of the proteins which may be reﬁned by sequencing RNAs from other tissues. secreted in saliva and their gene expression level in SGs Nonetheless, salivary genes not up-regulated in SGs should may result from insufﬁcient statistical power and/or regula- not be ignored. Some may indeed encode salivary effectors tion of salivary protein secretion. Plant cues may be required like Mif1, which suppresses induction of plant defense to trigger certain aphid responses, including protein Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 1567 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Boulain et al. GBE secretion (Powell et al. 2006), and this may explain why SG-expressed ones, especially when considering that some proteins encoded by highly transcribed salivary genes 90% of the proteins detected in saliva in this study were were not detected in the artiﬁcial diet. For example, C002 is encoded by SG-up genes. These SG-up effector candidates reported to be secreted into the plant and required for are more likely secreted by the aphid and thus exposed to aphid feeding (Mutti et al. 2008), but was not detected selective pressures exerted by the plant surveillance and by our saliva proteomics despite its high transcription level. defense systems. Moreover, we did ﬁnd traces of positive The comparison of transcriptomic and proteomic results selection acting on certain sites and branches of SG-up revealed limits in the secretion prediction pipeline as 14 candidate effector genes. Based on these observations proteins detected in saliva were not predicted to be se- and previous functional studies, we argue that the evolu- creted. Five of them may be miss-annotated, but the rest tionary history of some salivary effectors has been might have been secreted by pathways that were not con- conditioned by the specialization of SG tissue during es- sidered for the secretion prediction. tablishment of aphid-plant interactions. Fast Evolutionary Rates of Candidate Effectors Genes Salivary Gene Family Expansions SG-up genes showed lower degree of gene conservation In addition to single-copy genes, evolutionary histories of than SG-expressed genes. Indeed, 50% of the A. pisum three multigenic families were examined. Remarkably, these SG-up effector set had either no ortholog in other genomes three gene families encoded nearly half of the proteins or only in Aphidomorpha and Aphididae, whereas 64.2% detected in saliva and showed gene expansion in the A. pisum of the SG-expressed effector set had orthologs found lineage. Irrespective of their speciﬁcity to aphid lineages (CRP among arthropods. Furthermore, single-copy SG-up genes family) or conservation among insects (ACE and apN), these presented higher dN/dS ratio than SG-expressed genes on families show the highest number of members in A. pisum, average. This indicates that the evolution of SG-up genes with 15, 10, and 47 genes, respectively. Gene duplication and tends to be less constrained (relaxed selection) or acceler- diversiﬁcation therefore appear to have largely contributed ated (positive selection). We did not observe such contrast to the battery of A. pisum salivary proteins. Accordingly, SG- in between AT-up and AT-expressed gene sets, highlighting up and SG-expressed effector sets were enriched in dupli- the peculiar evolutionary history of genes constituting the cated genes. Aphididae-speciﬁc duplications were also pre- SG-up effector set. sent, but in lower numbers, in the D. noxia genome for CRP Interestingly, the analysis of the functionally characterized and apN families, as well as in A. glycines for the CRP. single-copy genes revealed different selective regimes. The Interestingly, clade 4 of the apN family, which contains a more conserved effectors (e.g., Mif1, Armet,and Mp10) majority of candidate A. pisum SG-up effector genes, did with low dN/dS ratio (0.12) could be involved in fundamen- not comprise any A. glycines members. Further analyses tal functions required for plant feeding, for example, preven- may indicate whether the loss happened in an ancestor of tion of phloem clogging, repression/manipulation of general the Aphidini subtribe or only in A. glycines (clade 4 is present plant defenses, such that orthologous proteins of similar in Aphidomorpha) or resulted from technical artifacts. sequences may be effective on numerous host species (Bos Positive selection was detected on speciﬁc duplicated copies et al. 2010; Furch et al. 2015; Naessens et al. 2015). In con- in D. noxia, but a tissue speciﬁc expression analysis or saliva trast, 11 effectors with dN/dS ratios above 0.40 have been proteomics is required to check whether they are SG spe- characterized as modulators of plant-aphid interactions in cialized potential effectors. functional studies, and their non-synonymous divergence Effector genes are often aggregated in large clusters in the may reﬂect adaptation to speciﬁc host plants. Some of them genome of plant pathogens (e.g., CRN and RXLR effectors of (C002, Mp1 and Mp2) were indeed shown to act in a species- the oomycete Phytophthora infestans, Haas et al. 2009). speciﬁc manner to promote aphid performances on host These clusters are thought to result from non-allelic homolo- plants (Pitino and Hogenhout 2013). Moreover, sites detected gous recombination and tandem duplications (facilitated in as under positive selection in Mp1 (Pitino and Hogenhout repeat-rich genomic regions), generally associated with rapid 2013; this study) and Me10 salivary effector genes may be birth and death evolution (Jiang et al. 2008; Haas et al. 2009). involved in adaptation of their respective aphid lineages. Effector gene clusters have also been observed in a few her- Although we cannot globally evaluate the relative con- bivore insects, including the gall midge Mayetiola destructor, tributions of relaxed and positive selection regimes on which shows a massive expansion of effectors organized in faster gene evolution, several evidences underline the im- clusters (Zhao et al. 2015). Although some members of each portance of positive selection in shaping the evolution of of the three A. pisum gene families are clustered on the same some candidate effectors identiﬁed in our study. It would genomic scaffolds (supplementary table S6, Supplementary be unintuitive to hypothesize that SG-up candidate effec- Material online), we were not able to identify clear gene clus- tor genes are more prone to relaxed selection than ters resulting from tandem duplication events, possibly due to 1568 Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Evolution and Gene Family Expansions of Aphid Salivary Effectors GBE assembly errors in the pea aphid reference genome (Jaquiery of salivary gene families reported here in A. pisum and sub- et al. 2018). sequent putative positive selection on some gene copies high- The functions of the CRP, ACE and apN gene families in light the importance of gene duplication in adaptation in the insects, including aphids, are largely unknown. The A. pisum ancestry. Notably, A. pisum constitutes a complex Aphididae-speciﬁc CRP gene family was ﬁrst described by that comes from the recent adaptive radiation of least 15 Guo et al. (2014), who demonstrated that CRP-13 expression biotypes (races or cryptic species), each specialized to one or in A. pisum was induced by feeding on plant in comparison to few related plant species within the Fabaceae (legume) family feeding on an artiﬁcial diet. However, silencing of CRP-13 did (Peccoud et al. 2009a, 2009b; Nouhaud et al. 2014; Peccoud not affect aphid survival on host plant. Proteins with numer- et al. 2015). Population genomics analyses on A. pisum bio- ous cysteine residues, like those of the CRP family, were pre- types have highlighted salivary and chemosensory genes as viously characterized for their antifungal activities in different likely contributing to host-speciﬁc adaptation in the A. pisum insect species (Banzet et al. 2002). Interestingly, some small complex (Jaquiery et al. 2012; Smadja et al. 2012; Nouhaud secreted cysteine-rich proteins act as effectors in the interac- et al. 2014; Duvaux et al. 2015; Eyres et al. 2017). We may tion between the Asian soybean rust fungus, Phakopsora relate these results to the existence of many SG-up effector pachyrhizi, and its host plant (Qi et al. 2016). One of them candidates that have undergone duplications and episodic was shown to suppress plant immunity by interacting with a positive selection in the A. pisum lineage. A large repertoire soybean transcription factor essential for negative regulation of effector genes might have increased the chance of adap- of immunity. Thus, the various cysteine-rich proteins appear tive mutations selected during biotype formation, reﬂecting to be active in different systems including plant-pathogen evolutionary patterns seen in other parasites’ effectors interactions, so the aphid-speciﬁc CRP family may play a (Stergiopoulos et al. 2012; Zhao et al. 2015). role in aphid nutrition once injected into the host plant. Strikingly, we found only one member of each CRP and The A. pisum ACE1 and ACE2 genes were previously apN (clade 4) family in the genome of M. persicae, consistent reported to contribute to aphid growth on plants, but not with the overall lower rate of gene expansion compared with on artiﬁcial diet (Wang et al. 2015b). As the encoded proteins the A. pisum lineage (Mathers et al. 2017). M. persicae is are predicted to have protease functions, their involvement in perhaps the most generalist aphid known, being able to col- cleavage of certain plant defense proteins or signaling com- onize hundreds of plant species across 40 families (Blackman ponents is speculated (Wang et al. 2015b). Aminopeptidases and Eastop 2000). Monitoring of gene expression during host have been found in aphid guts and reported as involved in switches of M. persicae laboratory lineages have suggested digestion (Rahbe et al. 1995; Cristofoletti et al. 2003)or in that acclimation to different plant species was enabled by a virus binding (Linz et al. 2015), but little is known about their rapid transcriptional plasticity of duplicated genes (Mathers function in aphid saliva. Furch et al. (2015) showed that A. et al. 2017). By contrast, very little changes in gene expression pisum saliva was able to degrade a major phloem protein 1 have been measured in A. pisum lineages when shifted from (PP1) in vitro. Because PP1 is involved in protein deposition on one suitable legume plant to another (Eyres et al. 2016). sieve plates after severe metabolic disturbance (Gaupels et al. Therefore, while large gene families, including salivary effec- 2008), its degradation by aphid saliva suggests that proteases tors, may be key in the rapid diversiﬁcation of specialized A. in aphid saliva degrade PP1 to prevent sieve-element occlu- pisum on various plant species, transcriptional plasticity may sion triggered by aphid feeding (Furch et al. 2015). In all cases, enable ecological generalism in M. persicae. These hypotheses protease function relies on peptide processing activity ensured require further investigation. In particular, comparative geno- by the M2 peptidase domain. Despite the detection of positive mics using more species of the M. persicae and A. pisum selection in several A. pisum M1 and M2 peptidase genes, the lineages in combination with functional analyses of candidate active site HEXXH is conserved in most of the members, indi- effectors, would help to precisely date the origins of genetic cating a functional cleaving activity. Some copies may have lost changes and to establish more robust correlations between active sites or domains, particularly in the A. pisum apN family these changes and ecological shifts. members that were too divergent to be clustered, reﬂecting possible functional changes after high diversiﬁcation. Supplementary Material Gene family expansion was previously described in A. Supplementary data are available at Genome Biology and pisum, particularly in cathepsin B gut proteases (Rispe et al. Evolution online. 2007), chemoreceptors (Smadja et al. 2009), small RNAs ma- chinery (Jaubert-Possamai et al. 2010), amino acid transport- ers (Price et al. 2011) and cuticular proteins (The International Acknowledgments Aphid Genomics Consortium 2010). In some cases, fast evo- lution and positive selection occurred on recently duplicated This work is part of H.B. PhD thesis funded by Plant Health genes (Rispe et al. 2007; Smadja et al. 2009). The expansion and Environment division of INRA and Region Bretagne. This Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 1569 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Boulain et al. GBE Conant GC, Wolfe KH. 2008. Turning a hobby into a job: how duplicated work was funded by Agence Nationale de la Recherche (ANR) genes ﬁnd new functions. Nat Rev Genet. 9(12):938–950. Bugspit (ANR-13-JSV7-0012-01) to A.S. and ANR Speciaphid Cox J, et al. 2011. Andromeda: a peptide search engine integrated into the (ANR-11-BSV7-005-01) to J.C.S., Campus France Ulysses MaxQuant environment. J Proteome Res. 10(4):1794–1805. (907424D) to A.S. and J.C. and IOS-1258028 from the Cristofoletti PT, Ribeiro AF, Deraison C, Rahbe Y, Terra WR. 2003. Midgut National Science Foundation and Kansas State University to adaptation and digestive enzyme distribution in a phloem feeding insect, the pea aphid Acyrthosiphon pisum. J Insect Physiol. F.W. and J.O. The Q-Exactive quantitative mass spectrometer 49(1):11–24. was funded under the SFI Research Infrastructure Call 2012, Dobin A, et al. 2013. STAR: ultrafast universal RNA-seq aligner. Grant Number 12/RI/2346 (3) to Sean Doyle. We are grateful Bioinformatics 29(1):15–21. to the Phylloxera Genomic Initiative (and especially Franc¸ois Duvaux L, et al. 2015. Dynamics of copy number variation in host races of Delmotte, Claude Rispe and Denis Tagu) for providing the pea aphid. Mol Biol Evol. 32(1):63–80. Dyrløv Bendtsen J, Nielsen H, von Heijne G, Brunak S. 2004. Improved Daktulosphaira vitifoliae genome sequences. Funding for prediction of signal peptides: signalP 3.0. J Mol Biol. 340(4):783–795. Daktulosphaira vitifoliae clone Pcf genomic sequencing was Elzinga DA, Jander G. 2013. The role of protein effectors in plant–aphid provided by INRA (AIP Bioressources) and BGI Biotech in the interactions. Curr opin Plant Biol. 16(4):451–456. frame of i5k initiative. Parts of the D. vitifoliae transcriptomic Elzinga DA, De Vos M, Jander G. 2014. Suppression of plant defenses by a resources were obtained within the 1KITE projects (Bernhard Myzus persicae (green peach aphid) salivary effector protein. Mol Plant Microbe Interact. 27(7):747–756. Misof, Bonn, Germany). We thank John Reese (Kansas State Eyres I, et al. 2017. Targeted re-sequencing conﬁrms the importance of University) for providing R. padi and S. graminum aphid sam- chemosensory genes in aphid host race differentiation. Mol Ecol. ples. We thank Gaetan Denis, Jean-Franc¸ois Le Gallic, 26(1):43–58. Frederique Maheo, and Sylvie Tanguy for technical support. Eyres I, et al. 2016. Differential gene expression according to race and host We are grateful to Claude Rispe for critical reading of the M.S. plant in the pea aphid. Mol. Ecol 25(17):4197–4215. Febvay G, Delobel B, Rahb e Y. 1988. Inﬂuence of the amino acid balance on the improvement of an artiﬁcial diet for a biotype of Acyrthosiphon Literature Cited pisum (Homoptera: Aphididae). Can J Zool. 66(11):2449–2453. Forister ML, et al. 2015. The global distribution of diet breadth in insect Abdellatef E, et al. 2015. Silencing the expression of the salivary sheath herbivores. Proc Natl Acad Sci U S A. 112(2):442–447. protein causes transgenerational feeding suppression in the aphid Furch ACU, van Bel AJ, Will T. 2015. Aphid salivary proteases are capable Sitobion avenae. Plant Biotechnol J. 13(6):849–857. of degrading sieve-tube proteins. J Exp Bot. 66(2):533–539. Anisimova M, Yang Z. 2007. Multiple hypothesis testing to detect lineages Gaupels F, et al. 2008. Nitric oxide generation in Vicia faba phloem cells under positive selection that affects only a few sites. Mol Biol Evol. reveals them to be sensitive detectors as well as possible systemic 24(5):1219–1228. transducers of stress signals. New Phytol. 178(3):634–646. Atamian HS, et al. 2013. In planta expression or delivery of potato aphid Go ¨ tz S, et al. 2008. High-throughput functional annotation and data min- Macrosiphum euphorbiae effectors Me10 and Me23 enhances aphid ing with the Blast2GO suite. Nucleic Acids Res. 36(10):3420–3435. fecundity. Mol Plant Microbe Interact. 26(1):67–74. Guo K, et al. 2014. Characterization of an aphid-speciﬁc, cysteine-rich Banzet N, et al. 2002. Expression of insect cystein-rich antifungal peptides protein enriched in salivary glands. Biophys Chem. 189:25–32. in transgenic tobacco enhances resistance to a fungal disease. Plant Guy E, et al. 2016. Optimization of agroinﬁltration in Pisum sativum pro- Sci. 162(6):995–1006. vides a new tool for studying the salivary protein functions in the pea Bendtsen JD, Jensen LJ, Blom N, von Heijne G, Brunak S. 2004. Feature- aphid complex. Front Plant Sci. 7:1171. based prediction of non-classical and leaderless protein secretion. Haas BJ, et al. 2009. Genome sequence and analysis of the Irish potato Protein Eng Des Sel. 17(4):349–356. famine pathogen Phytophthora infestans. Nature 461(7262): Bent AF, Mackey D. 2007. Elicitors, effectors, and R genes: the new par- 393–398. adigm and a lifetime supply of questions. Annu Rev Phytopathol. Harris KF, Maramorosch K. 2014. Aphids as virus vectors. Amsterdam: 45:399–436. Elsevier. Blackman RL, Eastop VF. 2000. Aphids on the world’s crops: an identiﬁ- Hubner NC, et al. 2010. Quantitative proteomics combined with BAC cation and information guide. United Kingdom: John Wiley & Sons TransgeneOmics reveals in vivo protein interactions. J Cell Biol. Ltd. 189(4):739–754. Bos JIB, et al. 2010. A functional genomics approach identiﬁes candidate Janz N. 2011. Ehrlich and Raven Revisited: mechanisms underlying codi- effectors from the aphid species Myzus persicae (green peach aphid). versiﬁcation of plants and enemies. Annu Rev Ecol Evol Syst. PLoS Genet. 6(11):e1001216. 42(1):71–89. Camacho C, et al. 2009. BLASTþ: architecture and applications. BMC Jaqui ery J, et al. 2018. Disentangling the causes for faster-X evolution in Bioinformatics 10:421. aphids. Genome Biol Evol. 10:507–520. Carolan JC, et al. 2011. Predicted effector molecules in the salivary secre- Jaqui ery J, et al. 2012. Genome scans reveal candidate regions involved in tome of the pea aphid (Acyrthosiphon pisum): a dual transcriptomic/ the adaptation to host plant in the pea aphid complex. Mol Ecol. proteomic approach. J Proteome Res. 10(4):1505–1518. 21(21):5251–5264. Carolan JC, Fitzroy CIJ, Ashton PD, Douglas AE, Wilkinson TL. 2009. Jaubert-Possamai S, et al. 2010. Expansion of the miRNA pathway in the The secreted salivary proteome of the pea aphid Acyrthosiphon hemipteran insect Acyrthosiphon pisum. Mol Biol Evol. pisum characterised by mass spectrometry. Proteomics 9(9): 27(5):979–987. 2457–2467. Jiang RHY, Tripathy S, Govers F, Tyler BM. 2008. RXLR effector reservoir in Chaudhary R, Atamian HS, Shen Z, Briggs SP, Kaloshian I. 2014. GroEL two Phytophthora species is dominated by a single rapidly evolving from the endosymbiont Buchnera aphidicola betrays the aphid by superfamily with more than 700 members. Proc Natl Acad Sci U S A. triggering plant defense. Proc Natl Acad Sci U S A. 111(24): 105(12):4874–4879. 8919–8924. 1570 Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Evolution and Gene Family Expansions of Aphid Salivary Effectors GBE Jones P, et al. 2014. InterProScan 5: genome-scale protein function clas- Peccoud J, Simon J-C, McLaughlin HJ, Moran NA. 2009. Post-Pleistocene siﬁcation. Bioinformatics 30(9):1236–1240. radiation of the pea aphid complex revealed by rapidly evolving endo- Katoh K, Kuma K, Toh H, Miyata T. 2005. MAFFT version 5: improvement symbionts. Proc Natl Acad Sci U S A. 106(38):16315–16320. in accuracy of multiple sequence alignment. Nucleic Acids Res. Petersen TN, Brunak S, von Heijne G, Nielsen H. 2011. SignalP 4.0: dis- 33(2):511–518. criminating signal peptides from transmembrane regions. Nat Kondrashov FA. 2012. Gene duplication as a mechanism of genomic ad- Methods 8(10):785. aptation to a changing environment. Proc Biol Sci Pierleoni A, Martelli PL, Casadio R. 2008. PredGPI: a GPI-anchor predictor. 279(1749):5048–5057. BMC Bioinformatics 9:392. Kriventseva EV, et al. 2015. OrthoDB v8: update of the hierarchical catalog Pitino M, Coleman AD, Maffei ME, Ridout CJ, Hogenhout SA. 2011. of orthologs and the underlying free software. Nucleic Acids Res. Silencing of aphid genes by dsRNA feeding from plants. PLoS One 43(D1):D250–D256. 6(10):e25709. Krogh A, Larsson B, von Heijne G, Sonnhammer ELL. 2001. Predicting Pitino M, Hogenhout SA. 2013. Aphid protein effectors promote aphid transmembrane protein topology with a hidden markov model: appli- colonization in a plant species-speciﬁc manner. Mol Plant Microbe cation to complete genomes. J Mol Biol. 305(3):567–580. Interact. 26(1):130–139. Law CW, Alhamdoosh M, Su S, Smyth GK, Ritchie ME. 2016. RNA-seq Pohlert T. 2014. The Pairwise Multiple Comparison of Mean Ranks analysis is easy as 1-2-3 with limma, Glimma and edgeR. Package (PMCMR). R package. http://CRAN.R-project.org/package= F1000Research 5:1408. PMCMR, last accessed May 22, 2018. Law CW, Chen Y, Shi W, Smyth GK. 2014. voom: precision weights un- Powell G, Tosh CR, Hardie J. 2006. Host plant selection by aphids: behav- lock linear model analysis tools for RNA-seq read counts. Genome Biol. ioral, evolutionary, and applied perspectives. Annu Rev Entomol. 15(2):R29. 51:309–330. Letunic I, Bork P. 2007. Interactive Tree Of Life (iTOL): an online tool for phy- Price DRG, Duncan RP, Shigenobu S, Wilson ACC. 2011. Genome expan- logenetic tree display and annotation. Bioinformatics 23(1):127–128. sion and differential expression of amino acid transporters at the Liao Y, Smyth GK, Shi W. 2014. featureCounts: an efﬁcient general pur- aphid/Buchnera symbiotic interface. Mol Biol Evol. 28(11):3113–3126. pose program for assigning sequence reads to genomic features. Qi M, et al. 2016. A small cysteine-rich protein from the Asian soybean rust Bioinformatics 30(7):923–930. fungus, Phakopsora pachyrhizi, suppresses plant immunity. PLoS Linz LB, Liu S, Chougule NP, Bonning BC. 2015. In vitro evidence supports Pathog. 12(9):e1005827. membrane alanyl aminopeptidase N as a receptor for a plant virus in Rahb e Y, Sauvion N, Febvay G, Peumans WJ, Gatehouse AMR. 1995. the pea aphid vector. J Virol. 89(22):11203–11212. Toxicity of lectins and processing of ingested proteins in the pea aphid Lo ¨ ytynoja A, Goldman N. 2005. An algorithm for progressive multiple Acyrthosiphon pisum. Entomol Exp Appl. 76(2):143–155. alignment of sequences with insertions. Proc Natl Acad Sci U S A. R-Core Team. 2017. R: A Language and Environment for Statistical 102(30):10557–10562. Computing. https://www.R-project.org/, last accessed June 12, 2018. Mathers TC, et al. 2017. Rapid transcriptional plasticity of duplicated gene Rispe C, et al. 2007. Large gene family expansion and variable selective clusters enables a clonally reproducing aphid to colonise diverse plant pressures for cathepsin B in aphids. Mol Biol Evol. 25(1):5–17. species. Genome Biol. 18:27 Ritchie ME, et al. 2015. limma powers differential expression analyses for Moreno A, et al. 2011. Aphids secrete watery saliva into plant tissues from RNA-sequencing and microarray studies. Nucleic Acids Res. the onset of stylet penetration. Entomol Exp Appl. 139(2):145–153. 43(7):e47–e47. Mullen SP, Shaw KL. 2014. Insect speciation rules: unifying concepts in Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor speciation research. Annu Rev Entomol. 59(1):339–361. package for differential expression analysis of digital gene expression Mutti NS, et al. 2008. A protein from the salivary glands of the pea aphid, data. Bioinformatics 26(1):139–140. Acyrthosiphon pisum, is essential in feeding on a host plant. Proc Natl Robinson MD, Oshlack A. 2010. A scaling normalization method for dif- Acad Sci U S A. 105(29):9965–9969. ferential expression analysis of RNA-seq data. Genome Biol. 11(3):R25. Mutti NS, Park Y, Reese JC, Reeck GR. 2006. RNAi Knockdown of a salivary Rodriguez PA, Escudero-Martinez C, Bos JIB. 2017. an aphid effector transcript leading to lethality in the pea aphid, Acyrthosiphon pisum.J targets trafﬁcking protein VPS52 in a host-speciﬁc manner to promote Insect Sci. 6:1–7. virulence. Plant Physiol. 173(3):1892–1903. Naessens E, et al. 2015. A secreted MIF cytokine enables aphid feeding Sela I, Ashkenazy H, Katoh K, Pupko T. 2015. GUIDANCE2: accurate de- and represses plant immune responses. Curr Biol. 25(14):1898–1903. tection of unreliable alignment regions accounting for the uncertainty Nicholson SJ, et al. 2015. The genome of Diuraphis noxia, a global aphid of multiple parameters. Nucleic Acids Res. 43(W1):W7–W14. pest of small grains. BMC Genomics 16:429. Smadja C, Shi P, Butlin RK, Robertson HM. 2009. Large gene family Nouhaud P, et al. 2014. Genomic regions repeatedly involved in diver- expansions and adaptive evolution for odorant and gustatory recep- gence among plant-specialized pea aphid biotypes. J Evol Biol. tors in the pea aphid, Acyrthosiphon pisum. Mol Biol Evol. 27(9):2013–2020. 26(9):2073–2086. Pan Y, Zhu J, Luo L, Kang L, Cui F. 2015. High expression of a unique aphid Smadja CM, et al. 2012. Large-scale candidate gene scan reveals the role protein in the salivary glands of Acyrthosiphon pisum. Physiol Mol of chemoreceptor genes in host plant specialization and speciation in Plant Pathol. 92:175–180. the pea aphid. Evolution 66(9):2723–2738. Peccoud J, Mah eo F, de la Huerta M, Laurence C, Simon J-C. 2015. Stamatakis A. 2006. RAxML-VI-HPC: maximum likelihood-based phyloge- Genetic characterisation of new host-specialised biotypes and novel netic analyses with thousands of taxa and mixed models. associations with bacterial symbionts in the pea aphid complex. Insect Bioinformatics 22(21):2688–2690. Conserv Divers. 8(5):484–492. Stergiopoulos I, et al. 2012. In silico characterization and molecular evo- Peccoud J, et al. 2010. Evolutionary history of aphid-plant associations and lutionary analysis of a novel superfamily of fungal effector proteins. their role in aphid diversiﬁcation. C R Biol. 333(6–7):474–487. Mol Biol Evol. 29(11):3371–3384., Peccoud J, Ollivier A, Plantegenest M, Simon J-C. 2009. A continuum of The International Aphid Genomics Consortium 2010. Genome genetic divergence from sympatric host races to species in the pea sequence of the pea aphid Acyrthosiphon pisum. PLoS Biol. aphid complex. Proc Natl Acad Sci U S A. 106(18):7495–7500. 8(2):e1000313. Genome Biol. Evol. 10(6):1554–1572 doi:10.1093/gbe/evy097 Advance Access publication May 18, 2018 1571 Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018 Thorpe P, Cock PJ, Bos J. 2016. Comparative transcriptomics and proteo- Yang Z. 2007. PAML 4: phylogenetic analysis by maximum likelihood. Mol mics of three different aphid species identiﬁes core and diverse effec- Biol Evol. 24(8):1586–1591. tor sets. BMC Genomics 17(1):172. Yang Z, Nielsen R. 2000. Estimating synonymous and nonsynonymous Varden FA, De la Concepcion JC, Maidment JH, Banﬁeld MJ. 2017. substitution rates under realistic evolutionary models. Mol Biol Evol. Taking the stage: effectors in the spotlight. Curr Opin Plant Biol. 17(1):32–43. 38:25–33. Yang Z, Nielsen R. 2002. Codon-substitution models for detecting molec- Wang W, et al. 2015. Armet is an effector protein mediating aphid–plant ular adaptation at individual sites along speciﬁc lineages. Mol Biol Evol. interactions. FASEB J. 29(5):2032–2045. 19(6):908–917. Wang W, et al. 2015. Angiotensin-converting enzymes modulate aphid– Yang Z, Nielsen R, Goldman N, Pedersen A-MK. 2000. Codon-substitution plant interactions. Sci Rep. 5(1): 8885. models for heterogeneous selection pressure at amino acid sites. Wenger JA, et al. 2017. Whole genome sequence of the soybean aphid, Genetics 155:431–449. Aphis glycines. Insect Biochem Mol Biol. doi: 10.1016/j.ibmb.2017.01. Yang Z, Wong WSW, Nielsen R. 2005. Bayes Empirical Bayes inference of 005. amino acid sites under positive selection. Mol Biol Evol. Wiens JJ, Lapoint RT, Whiteman NK. 2015. Herbivory increases diversiﬁ- 22(4):1107–1118. cation across insect clades. Nat Commun. 6(1): 8370. Zhang Y, Fan J, Sun J, Chen J. 2015. Cloning and RNA interference analysis Will T, Vilcinskas A. 2015. The structural sheath protein of aphids is re- of the salivary protein C002 gene in Schizaphis graminum.J Integr quired for phloem feeding. Insect Biochem Mol Biol. 57:34–40. Agric. 14(4):698–705. Will T, Furch ACU, Zimmermann MR. 2013. How phloem-feeding insects Zhao C, et al. 2015. A massive expansion of effector genes underlies gall- face the challenge of phloem-located defenses. Front Plant Sci. 4:336. formation in the wheat pest Mayetiola destructor. Curr Biol. Wu J, Baldwin IT. 2010. New insights into plant responses to the attack 25(5):613–620. from insect herbivores. Annu Rev Genet. 44:1–24. Downloaded from https://academic.oup.com/gbe/article-abstract/10/6/1554/4999395 by Ed 'DeepDyve' Gillespie user on 20 June 2018

http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png

Genome Biology and Evolution Oxford University Press

http://www.deepdyve.com/lp/oxford-university-press/fast-evolution-and-lineage-specific-gene-family-expansions-of-aphid-KT9HrE4DU0

Fast Evolution and Lineage-Specific Gene Family Expansions of Aphid Salivary Effectors Driven by Interactions with Host-Plants

Loading next page...

Publisher: Oxford University Press
ISSN: 1759-6653
eISSN: 1759-6653
DOI: 10.1093/gbe/evy097
pmid: 29788052
Publisher site: See Article on Publisher Site

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Fast Evolution and Lineage-Specific Gene Family Expansions of Aphid Salivary Effectors Driven by Interactions with Host-Plants

Fast Evolution and Lineage-Specific Gene Family Expansions of Aphid Salivary Effectors Driven by Interactions with Host-Plants

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Fast Evolution and Lineage-Specific Gene Family Expansions of Aphid Salivary Effectors Driven by Interactions with Host-Plants

Fast Evolution and Lineage-Specific Gene Family Expansions of Aphid Salivary Effectors Driven by Interactions with Host-Plants

Abstract

Journal

Recommended Articles

There are no references for this article.

Our policy towards the use of cookies