{"paper_id":"4b592d81-6cdc-4d1a-9c19-0bb5650f38f8","body_text":"1 \n \nMulti-omics approaches define novel aphid effector candidates associated 1 \nwith virulence and avirulence phenotypes 2 \nPeter Thorpe1, Simone Altmann2, Rosa Lopez-Cobollo3, Nadine Douglas4,5, Javaid Iqbal3, Sadia Kanvil3, 3 \nJean-Christophe Simon6, James C. Carolan4, Jorunn Bos2,7*, Colin Turnbull3* 4 \n1Division of Computational Biology, School of Life Sciences, University of Dundee, Dundee, DD5 4EH , 5 \nUK; 2Division of Plant Sciences, S chool of Life Sciences, University of Dundee, Dundee, DD5 4EH, UK ; 6 \n3Department of Life Sciences, Imperial College London, London SW7 2AZ, UK; 4Department of Biology, 7 \nMaynooth University, Republic of Ireland ; 5School of Biology and Environmental Science, University 8 \nCollege Dublin, Dublin 2, Republic of Ireland;  6Institut de Génétique, Environnement et Protection des 9 \nPlantes (IGEPP), INRAE, 35653 Le Rheu, France; 7The James Hutton Institute, Invergowrie, Dundee DD2 10 \n5DA, UK. *Authors for correspondence: j.bos@dundee.ac.uk, c.turnbull@imperial.ac.uk.  11 \nABSTRACT 12 \nBackground. Compatibility between plant parasites and their hosts is genetically determined by both 13 \ninteracting organisms. For example, plants may carry resistance (R) genes or deploy chemical 14 \ndefences. Aphid saliva contains many proteins that are secreted into host tissues. Subsets of these 15 \nproteins are predicted to act as effectors, either subverting or triggering host immunity. However, 16 \nassociating particular effectors with virulence or avirulence outcomes presents challenges due to the 17 \ncombinatorial complexity. Here we use defined aphid and host genetics to test for co -segregation of 18 \nexpressed aphid transcripts and proteins with virulent or avirulent phenotypes. 19 \nResults. We compared virulent and avirulent pea aphid parental genotypes, and their bulk segregant 20 \nF1 progeny on Medicago truncatula genotypes carrying or lacking the RAP1 resistance quantitative 21 \ntrait locus. Differential gene expression analysis of whole body and head samples, in combination with 22 \nproteomics of saliva and salivary glands , enabled us to pinpoint proteins associated with 23 \nvirulence/avirulence phenotypes. There was relatively little impact of host genotype, whereas l arge 24 \nnumbers of transcripts and proteins were differentially expressed between parental aphids, likely a 25 \nreflection of their classification as divergent biotypes within the pea aphid species complex. Many 26 \nfewer transcripts intersected with the equivalent differential expression patterns in the bulked F1 27 \nprogeny, providing an effective filter for removing genomic background effects. Overall, there were 28 \nmore upregulated genes detected in the F1 avirulent dataset compared with the virulent one. Some 29 \ngenes were differentially expressed both in the transcriptome and in the proteome datasets, with 30 \naminopeptidase N prot eins being the most frequent  differentially expressed family. In addition, a 31 \nsubstantial proportion (27%) of salivary proteins lack annotations, suggesting that many novel 32 \nfunctions remain to be discovered.  33 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n2 \n \nConclusions. Especially when combined with tightly controlled genetics of both insect and host, multi-34 \nomics approaches are powerful tools for revealing and filtering candidate lists down to plausible genes 35 \nfor further functional analysis as putative aphid effectors. 36 \nKEYWORDS: Aphid, transcriptomics, proteomics, saliva, effector, virulence, avirulence 37 \n 38 \nBackground 39 \nCrop losses due to insect pests represent an enduring challenge for  agriculture and global food 40 \nsecurity. Aphids are a major problematic group, due both to the direct damage they cause by phloem 41 \nsap feeding and to indirect effects through acting as vectors for transmission of many viruses. Impacts 42 \nof pests are further exacerbated by the breakdown of genetically based crop resistance mechanisms 43 \ndue to selection pressures driving pest evolution, as well as evolved insecticide resistance. 44 \nIn contrast to related fields such as plant -pathogen interactions, the molecular relationships that 45 \ndetermine (in)compatibility of plant -aphid interactions are relatively poorly understood. Specific 46 \nresistance to plant pathogens frequently involves recognition of pathogen effectors, often by 47 \nresistance proteins (R) characterised by nucleotide -binding and leucine rich repeat (NLR) domains. 48 \nSeveral coiled coil domain NLR proteins have been implicated in resistance to aphids and their close 49 \nrelatives. For example, Mi -1, Vat and Bph14 confer resistance to certain biotypes of Macrosiphum 50 \neuphorbiae (potato aphid) [1], Aphis gossypii (melon-cotton aphid) [2] and Nilaparvata lugens (brown 51 \nplanthopper) [3], respectively. These NLR receptors are predicted to be involved in direct or indirect 52 \nrecognition of molecular signatures that insects, like plant pathogens, release inside their hosts. 53 \nIndeed, aphids secrete multiple effector  proteins into their saliva, that are then predicted to be  54 \ndelivered into plant tissues to modulate host cell processes and to suppress or trigger host defences 55 \n[4–7]. Although there is one recent report of the BISP effector from brown planthopper, an aphid 56 \nrelative, interacting with the BPH 14 NLR in rice [8], there are currently no examples where cognate 57 \naphid effector and NLR pairs have been fully defined. Improved molecular insights into virulence and 58 \nresistance mechanisms taking place during both compatible and incompatible plant-aphid interactions 59 \nare therefore a priority, and can provide essential knowledge for future development of durable aphid 60 \ncontrol strategies. 61 \nThe availability of extensive genome, transcriptome and resequencing resources for the model aphid 62 \nspecies Acyrthosiphon pisum (pea aphid) [9, 10]  have enabled comprehensive genome -wide 63 \nexplorations. There are also genomic sequences now available at NCBI and Aphid Base 64 \n(https://bipaa.genouest.org/is/aphidbase/) for more than 25 species of aphids and close relatives , 65 \noften associated with  gene predictions  and transcriptomes [11]. In addit ion, several papers have 66 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n3 \n \nattempted to define the aphid effectorome, either by direct analysis of salivary proteins, or by 67 \ntranscriptomics of salivary glands, coupled with filters for predicted secreted, non -trans-membrane 68 \nproteins [12–17]. Beyond the true aphids (superfamily Aphidoidea), there are now genomic resources 69 \nfor sister groups within the Hemiptera such as planthoppers, leafhoppers, psyllids, whitefly and scale 70 \ninsects (https://www.ncbi.nlm.nih.gov/assembly/?term=hemiptera) that likewise are major crop 71 \npests, alongside genomes for triatomines and bed bugs, hemipterans that feed on animal rather than 72 \nplant hosts. Outside the Hemiptera, genomic data have been published for sucking pests such as thrips 73 \nand spider mites that feed on plant tissues other than phloem  [18–20]. Genome, transcriptome and 74 \nproteome comparisons across clades may enable definition of putative effector subsets that are 75 \nnecessary for different feeding modes, and may provide insights into conserved and divergent modes 76 \nof action in terms of how the plant immune system is targeted to enable successful parasitism. 77 \nDespite the wide range of functional genomics studies published to date, one common limitation is 78 \nthe lack of understanding of the differences in effector complements between virulent (host -79 \ncompatible) and avirulent (host -incompatible) genotypes. Genetic differences operate at several 80 \ntaxonomic levels. First, there are major differences across aphid species in their host preferences and 81 \nhost compatibilities. Some species, such as peach potato aphid ( Myzus persicae) are generalists that 82 \ncan feed on at least 40 0 known plant species, making them widespread crop pests [21]. Others are 83 \nspecialists, such as pea aphid (A. pisum) that exclusively feeds on legumes (Fabaceae). Second, there 84 \nis substantial diversity within species such as A. pisum that has led to its description as a species 85 \ncomplex comprising several host races that each have a strong preference for particular legume 86 \nspecies, supported by robust molecular marker fingerprints for each host race [22, 23] . T here is 87 \nevidence of divergence and differential expression of chemosensory gene families such as odorant 88 \nreceptors across different pea aphid biotypes  [24, 25] , bu t causative relationships have yet to be 89 \nestablished for genes and proteins that govern the range of compatible and incompatible interactions 90 \nseen. There is also clear evidence that some host races can survive and sometimes thrive as migrants 91 \non hosts outside their preferred species range [22]. Finally, at the intra-specific level for both aphids 92 \nand hosts, there can be a wide range of compatibilities. For example, from testing eight genotypes of 93 \nA. pisum in combination with 23 different Medicago truncatula (Mt) accessions, we discovered high 94 \ndiversity in both species that did not correspond particularly strongly to host races or to geographic 95 \norigins of the host lines [26]. Parallel to this, crossing two divergent pea aphid biotypes to generate F1 96 \nrecombinant populations uncovered Mendelian segregation of virulence/avirulence on Mt genotypes 97 \ncarrying the RAP1 aphid resistance QTL [27, 28].  98 \nHere, we report global exploration of the molecular basis for aphid virulence and avirulence on 99 \ndefined host genotypes. Specifically, we aimed to link phenotypes to candidate effectors and related 100 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n4 \n \ngenes by multiple comparisons of the transcriptomes and proteomes of two divergent parental pea 101 \naphid clones, along with the transcriptomes of  segregating avirulent and virulent pool ed individuals 102 \nfrom within F1 cross populations (Fig. 1). We also critically analysed the effectiveness of combined 103 \nomics approaches as a means to robustly uncover proteins with pivotal biological roles, such as 104 \neffectors that determine the difference between virulent and avirulent outcomes.  105 \nResults and Discussion 106 \nGeneration and analysis of aphid populations for RNA-Seq analyses 107 \nIn our previous work [27], we had demonstrated Mendelian segregation of inheritance of virulent and 108 \navirulent phenotypes in F1 pea aphid populations derived from a cross between N116 and PS01 109 \n(virulent and avirulent parental clones, respectively) when infested on M. truncatula hosts carrying 110 \nthe RAP1 resistance QTL [28]. On this basis, we reasoned that the molecular basis of the difference 111 \nbetween virulent and avirulent aphids could be revealed by transcriptomic and proteomic analysis. 112 \nHowever, there were likely to be thousands of genetic and gene expression differences between the 113 \nparental genotypes, that are representatives of phenotypically contrasting biotypes within the highly 114 \ndiverse pea aphid species complex [22, 26] . This makes it difficult to discern unrelated gen omic 115 \nbackground differences from causative genes responsible for suppressing host immunity  or for 116 \ntriggering R -gene dependent defences. To address this  challenge, we  employed a bulk segregant 117 \nanalysis (BSA-) RNA-Seq approach that would both reduce the genetic background effects and allow 118 \nus to test for heritability of differentially expressed (DE) genes across parental and F1 generations. 119 \nEnabling this strategy first required us to re-create the segregating F1 populations previously reported 120 \n[27].  121 \nWe induced sexual forms of PS01 and N116 and conducted reciprocal crosses, leading to screening of 122 \na total of 78 F1 clones on two host plant genotypes carrying RAP1: Jemalong A17 (hereafter A17), the 123 \noriginal source of the identified RAP1 QTL, and a resistant near -isogenic line ( RNIL) derived from a 124 \nmapping population [29] using A17 as one of the parents. The RAP1 aphid resistance QTL is highly 125 \neffective against PS01 aphids , typically resulting in high mortality, whereas N116 aphids are 126 \nunaffected. Progeny were verified as true F1 hybrids by a panel of seven SSR markers [22] and by 127 \nscreening for maternal inheritance of secondary symbionts reported in the pea aphid  [30]. Using a 128 \nvirulence index based on a combination of aphid survival and reproduction, F1 clones were first ranked 129 \naccording to performance on A17 . Phenotypes ranged from fully virulent to fully avirulent 130 \n(Supplementary Material 1A), similar to previous findings [27], although in the present experiment the 131 \npopulation as a whole did not display complete segregation into discrete virulent and avirulent 132 \ncategories. As also previously shown, resistance in the RNIL was slightly weaker than in A17, with F1 133 \nclones ranging from virulent to avirulent, and importantly performance on the two host genotypes 134 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n5 \n \nwas significantly correlated (Pearson r 0.72, P 1.82 e-13). All F1 clones were virulent on hosts lacking 135 \nRAP1 (Supplementary Material 1B ). We then selected 22 sibling F1 clones from each end of the 136 \ndistribution to provide two bulk sample sets with the strongest virulent (VIR) and avirulent (AVR) 137 \nphenotypes for subsequent transcriptomic analysis. Fig. 2 shows the complete separation of the 138 \nselected clones into virulent and avirulent classifications. As a final check prior to RNA-Seq 139 \nexperiments, we re -confirmed separation of survival rates of these two subsets of clones on both 140 \nresistant host genotypes (Supplementary Material 1C). 141 \nTranscriptomic analyses 142 \nWe first ran a n RNA-Seq experiment using the parental clones N116 and PS01 infested onto either 143 \nA17 or the susceptible DZA315.16 host (hereafter DZA) for 24 h prior to collection of heads for RNA 144 \nextraction. The multiple aims were to enrich for transcripts from salivary glands that express candidate 145 \neffectors, to uncover the transcriptome differences between the parental aphid genotypes, and to 146 \nreveal the impact of host plant genotype. Each aphid x host combination was replicated three times, 147 \ngiving a total of 12 libraries, ranging from 6.8 to 10.6 million reads uniquely mapped to the reference 148 \ngenome (Supplementary Material 2A).  149 \nHierarchical c lustering and principal components analysis  (PCA) of the transcriptom ic expression 150 \nprofiles both indicated that the replicates of each treatment were closely correlated in all cases, so no 151 \ndatasets needed to be discarded (Fig. 3A,B). These analyses additionally revealed that samples were 152 \nseparated largely by aphid genotype rather than host plant treatment. Overall, the transcriptomes of 153 \nthe two aphid genotypes  on A 17 plants  were clearly differentiated , with a total of 483 genes 154 \nsignificantly upregulated in N116 and 452 in PS01 (log2 fold change >2.0, FDR <0.05; Supplementary 155 \nMaterial 3; Fig. 3C). Similarly, on DZA host plants, 395 and 363 genes were upregulated in N116 and 156 \nPS01, respectively. In contrast, expression of relatively few genes, between three and 27, across all 157 \nthe pairwise comparisons, was significantly affected by the host plant (Supplementary Material 3; Fig. 158 \n3C). Functions of the DE genes are considered below, in conjunction with the other transcriptomic and 159 \nproteomic experiments. 160 \nWe next undertook a larger RNA-Seq experiment, sampling whole aphid bodies in order to capture 161 \ntranscripts from all tissues. Using aphids infested onto A17 host plants for 24 h, w e again compared 162 \nN116 and PS01 parental clones, but this time alongside the bulked segregant pools of VIR and AVR F1 163 \nclones described above . Five biological replicates for each gave a total of 20 RNA libraries  each 164 \ncontaining 14 to 22 million reads that uniquely map to the reference genome (Supplementary Material 165 \n2B). 166 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n6 \n \nSimilar to  the heads experiment, multivariate analysis by hierarchical clustering and PCA both 167 \nindicated that all replicates within each sample type grouped together, and that each sample type was 168 \nclearly differentiated. As expected, the genetically divergent parents were again highly separated, 169 \nwhereas the two pooled F1 datasets were much closer to each other , as they contain 50% of each 170 \nparental genome, with each pool representing the average transcriptome of multiple independent F1 171 \nclones (Fig. 4A,B). 172 \nDifferentially expressed genes were identified for all pairwise comparisons between samples (Fig. 4C). 173 \nThe number of up and down-regulated genes between the parental pairs and the pair of F1 pools are 174 \nshown in Fig. 5A, with the gene lists provided in Supplementary Material 3 . Several hundred genes 175 \nwere differentially expressed in both the whole-body and head comparisons of the parents. Some of 176 \nthese DE genes likely reflect genomic differences between the parental clones that are representatives 177 \nof divergent pea aphid biotypes. However, relatively few DE genes were detected in the F1 samples, 178 \nwith only 24 genes up-regulated in the VIR pool and 64 in the AVR pool. These numbers can also be 179 \ninterpreted as a higher number of genes being down-regulated in the VIR F1 aphids. Fig. 5B,C show 180 \nthe overlaps across head and whole -body datasets for N116 /VIR and PS01 /AVR, respectively . 181 \nUnexpectedly, the intersections of DE genes revealed subsets where the direction of expression was 182 \nopposite between the parental pair and the F1 pooled pairs , with three genes upregulated in N116 183 \nand AVR F1, and 13 genes upregulated in PS01 and VIR F1 (Fig . 5D, Fig. 7G,H). Moreover, very few 184 \ngenes were upregulated in both parental N116 and VIR F1 pool datasets . A plausible explanation is  185 \nthat the genes governing virulence in N116 are not the same  as those that result in virulent 186 \nphenotypes in the F1 population . Each individual in the F1 population carries a random 50% of the 187 \ngenome of each parent , creating a high degree of combinatorial complexit y. Nonetheless, the DE 188 \ngenes in the F1 data derive from the average across the 22 individuals used to create each bulk RNA 189 \npool, and are therefore likely to be biologically relevant to virulence or avirulence functions rather 190 \nthan background genomic noise.  Such genes merit further exploration in both parental and F 1 191 \ngenotypes. 192 \nQuantitative proteomic analysis of saliva and salivary glands. 193 \nTo determine whether differences exist between the salivary protein profiles of the two parental 194 \naphid clones, a comparative analysis of salivary gland and salivary proteomes was conducted. A total 195 \nof 2343 and 2276 high confidence proteins were detected from salivary glands of N116 and PS01 , 196 \nrespectively (Supplementary Material 4 ), with 2038 proteins (80%) common to both ( Fig. 6A). Each 197 \nbiotype had similar proportions of non-annotated proteins (PS01: 5.4 % and N116: 6.2%) and proteins 198 \npredicted to have secretion signals (PS01: 16.6% and N116: 17.3%). These proportions of secreted and 199 \nnon-annotated proteins are typical for pea aphid biotypes [12, 31]. Two major clusters were revealed 200 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n7 \n \nby PCA (Fig. 6C), corresponding to the two aphid genotypes. Principal Components 1 and 2 account 201 \nfor 64% of the variation, indicating distinct protein profiles in the salivary glands of each genotype. 202 \nThis distinction was further supported by quantitative analysis that identified 23 5 statistically 203 \nsignificant differentially abundant (SSDA) proteins (p<0.05), with 1 36 and 99 proteins having higher 204 \nabundances in N116 and PS01 salivary glands , respectively ( Fig. 6E; Supplementary Material 4 ). 205 \nRelative fold changes (RFC) ranged from −48.5 to +140.0 indicating that even when both genotypes 206 \nengage in compatible interactions with the same plant type ( V. faba in this case) the salivary gland 207 \nprofiles are divergent both qualitatively and quantitatively.  208 \nOf the 136 SSDA salivary gland proteins with increased abundance in N116, 60 (44%) were predicted 209 \nto be secreted and 27 (20%) had no annotations. Similar proportions were observed within the 99 210 \nSSDA proteins with increased abundance in PS01 , with 33 (33%) and 18 (18%) proteins having a 211 \nsecretion signal or no annotations, respectively. The se proportions of secreted and non -annotated 212 \nproteins within the differentially abundant sets  are substantially higher than the corresponding 213 \nproportions in the background salivary gland proteomes described above. Of the top ten proteins with 214 \nthe highest relative abundance in N116, seven had no annotation : ACPISUM_000319 (ACYPI007553; 215 \nRFC 140.0) and ACPISUM_029783 ( LOC100573424; RFC 64), ACPISUM_008675 (LOC100162547; RFC 216 \n32), ACPISUM_016335 (Not annotated; RFC 26), ACPISUM_017388 (LOC103309964; RFC 21.1), 217 \nACPISUM_003551 (LOC100534636; RFC, 21.1) and ACPISUM_009099 (LOC112598674, 18.4). The 218 \nother proteins in the top ten were a kinase ACPISUM_015393 (developmentally -regulated protein 219 \nkinase 1; RFC 64) and two aminopeptidases (ACPISUM_009259; RFC 36.8 and ACPISUM_005699; RFC 220 \n22.6). Of the top ten proteins with highest abundances in PS01 in comparison to N116, two were 221 \nuncharacterised: ACPISUM_007394 (LOC100572241; RFC 48.5) and ACPISUM_007714 222 \n(LOC100534636; RFC 11.3) ; and two were glutathione S -transferases (ACPISUM_019160 and 223 \nACPISUM_001883, both RFCs of 8.6). Other proteins included a  different developmentally-regulated 224 \nprotein kinase (ACPISUM_005630; RFC 17.1), a peroxidase (ACPISUM_020816; RFC 9.8), a prostatic 225 \nspermine-binding protein (ACPISUM_004331; RFC 8), peroxidasin (ACPISUM_019870; RFC 6.5), an 226 \nATPase subunit (ACPISUM_009308; RFC 5.7) and glyoxylate reductase (ACPISUM_021751, RFC 4.9). 227 \nWe next examined aphid saliva proteins. Although the samples are collected from artificial diets, these 228 \nsalivary secretomes are likely to be highly similar to the proteins delivered into plant tissues during 229 \ninteractions with the host, and therefore are predicted to include the entire set of effectors. We 230 \nfocussed on categorisation of the total salivary protein lists, and of the DE proteins. Although the 231 \nanalysis of saliva revealed far fewer proteins than from the salivary gland samples, there is again a 232 \nclear distinction between the two genotypes. A total of 69 and 97 high confidence proteins were found 233 \nin N116 and PS01 saliva, respectively (Fig. 6B; Supplementary Material 4) with 22 (32% for N116) and 234 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n8 \n \n50 (52% for PS01) proteins, being deemed unique to each. A large proportion (30% for PS01 and 25% 235 \nfor N116) of the salivary proteomes had no annotations, indicating their potential phylogenetic 236 \nrestriction to aphids. In addition, 39% and 3 2% of the proteins had predicted canonical secretion 237 \nsignals for PS01 and N116 saliva, respectively. Notably, although s aliva proteins detected in diet 238 \nsamples have, by definition , been secreted , the majority appear not to have canonical secretion 239 \nsignals. Explanations range from incomplete/incorrect gene models to non-canonical or alternative 240 \nsecretion mechanisms. Our results highlight the importance of combining several approaches when 241 \nattempting to identify potential effectors and molecular determinants of virulence/avirulence. 242 \nOmitting proteins without secretion signals from bioinformatic pipelines may result in many effector 243 \ncandidates being overlooked.  244 \nAs with the salivary glands, PCA of the salivary proteins completely resolved two groups, with PC1 and 245 \nPC2 accounting for 94% of the total variation (Fig. 6D). Label free quantitative analysis using MaxQuant 246 \nidentified 47 SSDA proteins with 12 and 35 proteins having higher abundance in N116 and PS01 saliva, 247 \nrespectively ( Fig. 6F; Supplementary Material 4 ). Notably, N116 saliva comprises fewer detected 248 \nproteins and fewer SSDA proteins than PS01, possibly pointing to a strategy that enables evasion of 249 \nhost defences. If , for example, one or more of the proteins uniquely detected in PS01 saliva act as 250 \navirulence factors due to cognate receptors in the host plant, their absence or low abundance in N116 251 \nmay result in a compatible interaction. However, it remains to be experimentally determined whether 252 \nthese genotypic differences in type or number of saliva proteins are causatively associated with 253 \nvirulence or avirulence.  254 \nMost of the salivary proteins identified here have previously been associated with pea aphid saliva 255 \nincluding multiple members of M1 and M2 metalloprotease families , along with peroxidases, 256 \nglutathione-S-transferases, glucose dehydrogenase  and regucalcin [12, 32] . Apart from the 257 \nAminopeptidase N (APN) category discussed in detail below, the most frequent annotation was for 258 \nunknown proteins: 20-26% of the total  saliva list for each clone, and 21% of the DE saliva proteins. 259 \nFour out of the ten DE unknown proteins also featured within the top 20 proteins by MS intensity or 260 \nprotein coverage. High proportions of unknown proteins have been noted in earlier studies of aphid 261 \nsaliva and the salivary gland predicted secretome [31]. In addition, a homologue of a salivary effector 262 \npreviously characterised for Myzus persicae (Mp1) [33] had a higher abundance in PS01 saliva 263 \n(ACPISUM_000421; RFC 14). The relative fold changes of salivary proteins ranged from -2352 for 264 \nregucalcin to 724 for members of the APN (M1 metalloprotease) family, which represented the most 265 \ndifferentially abundant proteins in PS01 and N116 saliva, respectively. Although these RFC values can 266 \nbe considered arbitrary due to imputation of low abundant values in samples where the proteins are 267 \nin fact absent, the re is very clear divergence of salivary proteomes both in the proteins uniquely 268 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n9 \n \ndetected in one or other genotype , and in the large differences in apparent abundance of several 269 \nproteins present in both genotypes. The full lists of proteins exclusively found in the saliva or salivary 270 \ngland proteomes of both genotypes are provided in Suppleme ntary Material 4 , with 25 and five 271 \nproteins exclusive to the salivary glands and saliva of N116, respectively. For PS01, the corresponding 272 \nnumbers were  10 and 13 proteins exclusive to the salivary glands and saliva, respectively. These 273 \nproteins were present in all replicates of one genotype while being absent in all replicates of the other, 274 \nstrongly supporting their status as candidate effectors, that may individually or collectively determine 275 \nthe VIR and AVR phenotypes observed for each genotype on different host plants. 276 \nComparison of the quantitative differences in protein abundance across both the saliva and salivary 277 \ngland datasets revealed clear similarities in the two proteomes analysed for each genotype . Five  278 \nproteins that were of higher abundance in N116 saliva were also more abundant  in N116 salivary 279 \nglands in comparison to their PS01 counterparts. A similar trend was observed for nine PS01 salivary 280 \nand salivary gland proteins ( Supplementary Material 4 ), with the RFCs for these proteins positively 281 \ncorrelated across both biological sample types. The fact that the abundances of these salivary gland 282 \nproteins are mirrored at the level of externally delivered oral secretions highlights the  robustness of 283 \nboth analyses, and points to likely roles as virulence or avirulence determinants in two genotypes with 284 \ndistinct host preferences. Such proteins represent excellent candidates for future characterisation to 285 \ndetermine their effector status , especially those that are also supported by DE transcript profiles 286 \n(Table 1).  287 \nOverlap between transcriptomics and proteomics datasets  288 \nAcross the transcriptomics and prote omics experiments, we analysed all the intersections then 289 \nextracted the proteins and DE gene subsets that showed the greatest overlaps  (Table 1; 290 \nSupplementary Material 3 and 4), partitioning into genes/proteins associated with virulence, in N116 291 \nor the VIR F1 pool, or with avirulence, in PS01 or the AVR F1 pool. The number of DE genes or proteins 292 \nin the hea d transcriptome, whole body transcriptome and salivary gland proteome datasets were 293 \nbroadly similar between VIR and AVR samples. However, the PS01 saliva protein and the AVR F1 pool 294 \ntranscript lists were longer than those for N116 saliva and VIR F1 pool transcripts, reflected by larger 295 \nintersections in the former. Over half (33/64) of genes upregulated in the AVR F1 pool were also in at 296 \nleast one other list, whereas only three out of 24 intersected from the VIR F1 pool data. Whole body 297 \nRNA-Seq data for a selection of these intersected genes are plotted in Fig . 7. Several of the AVR -298 \nupregulated genes shown are annotated as enzymes with hydrolase, glycosidase or peroxidase 299 \nfunctions. Other annotations include a transcription factor and proteins of unknown function. Genes 300 \non the VIR side included  ACPISUM_013796 (myrosinase 1 -like) and ACPISUM_019971 (glutathione 301 \nhydrolase 1 proenzyme -like), although these were not found in saliva . Across the multiple 302 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n10 \n \nexperiments, the two most frequently found genes in the AVR data were ACPISUM_021997 303 \n(regucalcin-like) previously reported as a Ca-binding protein [32], present in all lists except heads RNA, 304 \nand ACPISUM_029930 (uncharacterized protein LOC100575698), present in all five lists. These AVR-305 \nrelated salivary proteins represent strong candidates for functional effectors, based on the multiple 306 \nstrands of evidence for their differential expression and importantly for co-segregation of their 307 \nexpression with the avirulence phenotype  in the F1 population . We have therefore uncovered 308 \nheritable differences in salivary proteins that associate with avirulence, in this case  an incompatible 309 \nphenotype on Mt hosts carrying the RAP1 QTL [27, 28]. Intriguingly, however, we found no equivalent 310 \nstrong candidates for salivary proteins that might represent the dominant virulence factor predicted 311 \nby previous genetic studies [27]. Alternative explanations for the Mendelian segregation found in that 312 \nstudy could be that the proposed “virulence” gene is not an effector per se, but instead  could be an 313 \nupstream positive regulator, or a negative regulator of one or more effectors that act as avirulence 314 \nfactors detected by a RAP1 dependent pathway.  315 \nGene Ontology analysis 316 \nWe undertook Gene Ontology (GO) analysis to reveal functional categories and genes that were 317 \nenriched in the differentially expressed gene and protein data sets. Using a FDR of <0.05, many gene 318 \nsets contained few or no significant ly enriched terms (Table 2; Supplementary Material 5). For the 319 \nwhole-body transcriptome data, aminopeptidase N (APN) proteins were strongly enriched, with 320 \ndifferent genes within this family upregulated in each of the parental aphids (discussed further below). 321 \nThese trends were reinforced by comparison of parental transcriptomes in the heads RNA-Seq 322 \nanalyses where APN proteins were similarly enriched in both parents. The DE gene sets between the 323 \npooled VIR and AVR F1 samples indicated no enriched terms in the VIR data, and only a single term 324 \namong the AVR upregulated genes: glucosidase II complex, that localises to the ER. These two gene 325 \nsets are both relatively small (64 and 24 genes), reducing the likelihood of finding significant trends. 326 \nBecause very few significantly enriched terms were revealed by the initial GO analyses, we applied a 327 \nlower stringency to inform wider trends in each of the DE gene sets. Here, we examined all terms for 328 \nwhich at least two genes and a significant P value (<0.05) were returned. For the DE gene sets from 329 \nRNA-Seq of heads, the majority of enriched terms were associated with the virulent N116 parent on 330 \nboth host genotypes. Although there was obvious redundancy of many terms, a substantial proportion 331 \n(30-40%) for N116 relate to energy metabolism including mitochondria, TCA cycle, oxidative 332 \nphosphorylation and lipid metabolism. In contrast, the PS01 enriched terms included several for 333 \nprotein processing including peptidases, proteolysis and protein glycosylation ; and several for ATP-334 \nrelated transport (Supplementary Material 5 ). When each parental aphid genotype was compared 335 \nseparately for its differential responses to the two host genotypes (A17 and DZA), no significant terms 336 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n11 \n \nwere found for PS01, and only one weakly significant term for N116: polytene chromosome puffing. 337 \nThe equivalent GO analysis of whole body RNA-Seq data returned significantly enriched terms for both 338 \naphid genotypes, including several for protein modification (Supplementary Material 5). 339 \nFor the DE datasets from salivary gland proteomes, the lower stringency analysis revealed enrichment 340 \nof distinct functional categories for each parental genotype. For N116, protein modification terms 341 \nwere prevalent including peptidase activity, serine -type endopeptidase inhibitor activity, negative 342 \nregulation of protein metabolic process, aminopeptidase activity, protein kinase binding and 343 \nregulation of protein phosphorylation. In contrast, for PS01, ATPase terms were predominant 344 \nincluding several related to membrane transport , as also found in the PS01 heads RNA -Seq data  345 \n(Supplementary Material 5).  346 \nExopeptidases are abundant in saliva, and the majority are DE between aphid genotypes 347 \nThe saliva protein total and DE lists were much shorter, precluding formal GO analysis, but manual 348 \ninspection indicated high proportions of exopeptidases: a total of 29 different proteins (Table 3) , 349 \nrepresenting 22-34% of the protein list for each genotype. These were mainly APN proteins but also  350 \nfour members of the Angiotensin Converting Enzyme (ACE) family that are M2 metalloproteases with 351 \ncarboxypeptidase activity. The abundance of APNs in the saliva protein list broadly corroborates the 352 \nmajor enriched GO categories detected in the transcriptome analyses. 353 \nMost of the exopeptidases detected from aphid saliva ( 23/29; 7 9%) were differentially abundant 354 \nbetween the parental aphid genotypes. Twenty-two of the 29 saliva exopeptidases were also found in 355 \nthe salivary gland proteomes, with many showing the same direction of differential expression (9 APN, 356 \n2 ACE). Moreover, 1 5 (60%) of th e APN proteins were DE in heads and /or whole body RNA-Seq 357 \nsamples (Table 3). Previous reports on pea aphid saliva and salivary gland components have also 358 \nreported multiple APN and ACE proteins [12, 13, 32, 34]. Similar to our findings, one of these studies 359 \nreported 11 APN genes that were differentially expressed in a biotype -specific manner, with five of 360 \nthese detected as proteins in saliva [13]. Taking all the evidence together, it is clear that the APN family 361 \nis highly diversified in pea aphids and represents a major component of the salivary proteome by 362 \nseveral measures: the high total number of proteins detected , many of these proteins are high 363 \nabundance (13 of 20 top scoring in both N116 and PS01 saliva), and most are differentially expressed 364 \nbetween aphid genotypes. 365 \nAphid and mammalian ACE proteins have similar sequences and may have broadly similar functions 366 \nas dipeptidases or by cleaving a single amino acid from the C terminus. However, mammalian ACE 367 \nproteins are membrane anchored whereas aphid ACEs carry secretion signals, consistent with their 368 \ndetection in saliva. The exact catalytic functions and biological roles of aphid ACE and APN proteins 369 \nremain to be determined. Cleavage of proteins and peptides could relate to targeting host proteins 370 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n12 \n \nsuch as those involved in defensive sieve -tube blocking as shown at least for the atypical 371 \nextrafascicular phloem exudate of cucurbits [35]. Alternatively, although there is currently no direct 372 \nevidence, exopeptidases may act on other salivary protein components, for example to process 373 \neffectors into active forms. Another non-mutually exclusive possibility is a role in aphid nutrition, with 374 \nmany insects using extra-organismal (extra-oral) digestion typical of arthropods including Hemiptera, 375 \nenabling nutrition capture from large host s/prey [32, 36] . Exopeptidases typically release N or C 376 \nterminal single amino acids and dipeptides, potentially enabling supply of essential amino acids, some 377 \nof which cannot be biosynthesised directly from the enzyme repertoires of hemimetabolous aphids.  378 \nMulti-omic approaches to detecting candidate effectors 379 \nWe compared the efficiencies of the four different experiments in terms of detecting aphid candidate 380 \neffectors and related genes: RNA-Seq of heads and whole bodies, and proteomics of saliva and salivary 381 \nglands. For all datasets, we focussed mainly on differential expression between the highly divergent 382 \nparental clones N116 and PS01. Because saliva represents the “ground truth” of proteins predicted to 383 \nbe delivered into plant host tissues, we additionally considered saliva proteins that were detected but 384 \nnot DE. Although the proteomics methods are highly sensitive, there are likely to be some further low 385 \nabundance salivary proteins that were not detected here. In addition, there may be some salivary 386 \nproteins that are only expressed in response to aphids interacting with their host plants , and hence 387 \nwould not be found in artificial diet samples. Likewise , some proteins may not be stable under the 388 \nartificial diet conditions. As a case study, we selected the significantly enriched exopeptidases , that 389 \ncomprised the large APN family and the smaller group of ACE proteins. We compared success of 390 \ndetecting genes from the saliva data in the other three experiments, and noted whether the same DE 391 \npatterns were found (Table 3). The overall trends were broadly correlated, with 18/24 (75%) DE saliva 392 \nproteins also found to be DE in at least one of the other approaches. Only two genes showed a 393 \nmismatch in DE direction: ACPISUM_009259 between salivary gland and whole body; and 394 \nACPISUM_020790 between saliva and salivary gland. Individually, RNA-Seq of heads was the most 395 \neffective experiment (14/24) at corroborating the DE saliva protein evidence, followed by RNA-Seq of 396 \nwhole bodies (10) and proteomics of salivary glands (8).  397 \nThere are several reports where effectors are predicted f rom aphid salivary gland transcriptomes or 398 \nproteomes, or other transcriptome datasets, typically filtering for presence of a signal peptide or other 399 \nsecretion motif, and absence of transmembrane domains [12–17]. For our exopeptidase data (Table 400 \n3), we detected an additional seven APNs in salivary gland proteomes or the transcriptome data, that 401 \nwere not found in saliva, of which five were DE in at least one dataset. The ir absence from saliva 402 \nindicates these proteins may be considered false positives for candidate effectors, although some low 403 \nexpressed proteins may go undetected . We considered which of the approaches was the most 404 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n13 \n \neffective at detecting candidate effectors, and whether multiple omics approaches are advantageous, 405 \nnoting that all require substantial resource investment. Although saliva collection is an exacting and 406 \ntime-intensive procedure, saliva proteomics provided the greatest coverage of candidate effectors 407 \nhere, and quantitative analysis of mass spectrometry data enables robust assignment of differential 408 \nexpression. Of the other approaches, RNA-Seq of heads may be the most effective means to 409 \ncomplement the saliva analyses by reinforcing evidence of differential expression, but in the work 410 \nhere did not greatly extend the effector lists per se.  411 \nConclusion 412 \nIn this study, we demonstrated that transcriptomics and proteomics are both highly effective tools for 413 \ndiscovering differentially expressed aphid genes and proteins . The protein subsets present in saliva 414 \nare likely candidates for effectors with virulence and/or avirulence functions in host plants , and 415 \nrepresent priorities for further study especially to determine if differential protein abundance is 416 \ninherited into the segregating F1 aphid populations. Precise biochemical functions and host targets of 417 \nmost of these effectors are also currently unknown even in cases, such as the exopeptidases, where 418 \nthere are confident gene annotations. Exopeptidases are dominant in saliva by number of different 419 \nproteins, by frequency of differential abundance, and by quantity. Likewise, there are many proteins 420 \nof unknown function, with a substantial proportion found at high  levels in saliva.  Some of these 421 \nunknown proteins may prove to be pivotal in explaining aphids’ unique and highly successful lifestyle 422 \nas phloem feeders. 423 \nMethods 424 \nAphids and crossing 425 \nPea aphid (Acyrthosiphon pisum) clones were maintained on tic bean (Vicia faba minor) as described 426 \nin [26]. Parental genotypes were PS01 and N116. PS01 is a biotype adapted to Pisum sativum whereas 427 \nN116 is adapted to Medicago sativa [26]. Reciprocal crosses were made between PS01 and N116 to 428 \ngenerate F1 hybrid populations, following the protocol of [27]. In brief, parthenogenetic females were 429 \ninduced to generate sexual forms by transfer to short days and lower temperatures to simulate 430 \nautumn. Eggs resulting from controlled matings were collected onto moist filter paper in petri dishes, 431 \nand subjected to 90 to 105 days at 4°C to induce exit from diapause. Individual hatchlings were 432 \nsubsequently used to generate multiple parallel clonal F1 lineages. Parents and progeny were 433 \ngenotyped with a set of seven microsatellite markers [22] to verify correctness of crosses. All new F1 434 \nprogeny were maintained for at least three generations before testing performance on different host 435 \nplants.  436 \nPlants and assessment of virulence 437 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n14 \n \nBased on previous findings [27], PS01 aphids are avirulent on Medicago truncatula J A17 that carries 438 \nthe resistance QTL, RAP1 [28]. Near isogenic lines (NILs) derived from a cross (LR4 [29]) between A17 439 \nand M. truncatula DZA315.16 were also used. PS01 is likewise incompatible with the resistant NIL 440 \n(RNIL), but is compatible with the susceptible NIL (SNIL) and with DZA315.16. N116 aphids are 441 \ncompatible with all these genotypes. F1 progeny were tested for virulence on both A17 and RNIL, 442 \nbased on [26]. Briefly, five nymphs of each clone were infested onto ten A17 or RNIL plants, then 443 \nscored for survival and production of new nymphs 10 d later. At least 40 F1 clones each of PS01 x N116 444 \nand N116 x PS01 were screened. An overall virulence index was adapted from a calculation proposed 445 \nin [37]: 446 \nVirulence index = log2 (mean number surviving out of 5 x number of nymphs produced + 1) 447 \nVirulent (VIR) clones were defined as index >4 and >5 on A17 and RNIL, respectively, and avirulent 448 \n(AVR) clones were correspondingly defined as index <2 and <4. The different category thresholds on 449 \nA17 and RNIL reflect the latter’s slightly lower resistance. Clones falling into the same phenotype 450 \ncategory (VIR or AVR) on both A17 and RNIL were then subject to a further confirmation screen where 451 \nsurvival on A17 and RNIL was counted 5 d after infestation. In the confirmation experiment, four plants 452 \nwere used for each aphid x host combination, with five aphids infested onto each plant. Cutoffs were 453 \n>80% survival for virulence on both hosts, and <20% and <70% for avirulence on A17 and RNIL, 454 \nrespectively. A few F1 clones showed relatively high survival at 5 days but had very weak growth, and 455 \ntherefore were categorised as AVR. Only F1 clones displaying the same phenotype category on all 456 \nscreening experiments were used subsequently in molecular experiments. 457 \nSampling for RNA-Seq 458 \nHeads experiment: Young adult aphids of clones N116 and PS01, cultured on Vicia faba minor, were 459 \ninfested onto either A17 or DZA315.16 M. truncatula plants for 24 h, then heads (40 per sample) were 460 \ndissected and frozen immediately on dry ice then stored at -80°C. Three replicates were done for each 461 \naphid x plant combination. 462 \nWhole body experiment: Samples were parental aphid clones (N116 and PS01) and pools of VIR and 463 \nAVR F1 progeny. Aphids of each individual genotype, age 2 to 3 d, were placed on independent A17 464 \nplants for 24 h then frozen in liquid nitrogen and stored at -80°C until processing. A total of 22 VIR and 465 \n22 AVR F1 aphid clones were collected individually, before pooling five aphids of each genotype  to 466 \ncomprise one  sample. Five biological replicates were analysed for both parental and pooled F 1 467 \ngenotypes. 468 \nRNA extraction, library preparation and sequencing 469 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n15 \n \nHeads were dissected and processed as described in [16]. Total RNA was extracted using a plant RNA 470 \nextraction kit (Sigma -Aldrich). Illumina TruSeq stranded mRNA -Seq libraries were sequenced at the 471 \nGenome Sequencing Unit at the University of Dundee on an Illumina HiSeq 2000. 472 \nRNA for the BSA-RNA-Seq analysis was isolated from three two to three day old nymphs of parental 473 \nlines (N116, PS01), 22 VIR F1 lines and 2 2 AVR F1 lines, using the Norgen Plant and Fungal RNA kit 474 \n(Sigma E4913). The RNA isolation followed the instructions of the company supplementing Lysis buffer 475 \nC with ß-mercaptoethanol. An on -column DNase digest was performed (RNase -Free DNase Set, 476 \nQiagen) and the concentration of each sample determined via a Qubit fluorometer with the QubitTM 477 \nRNA Broadrange (BR) assay kit (Thermo Fisher Scientific). Samples corresponding to five replicates of 478 \neach of the parental lines and the VIR and AVR F1 pools were used to generate a total of 20 Illumina 479 \nTruSeq stranded mRNA-Seq libraries which were sequenced in 150 bp paired-end mode on an Illumina 480 \nHiSeq4000 at Edinburgh Genomics. 481 \nRNA-Seq data processing and visualisation. 482 \nIllumina RNA sequence reads were subjected to quality control using FastQC. The reads were the 483 \ntrimmed using Trimmomatic (version 0.32) Q15, min length 55. The trimmed fastq files were the n 484 \nquasi mapped to the nucleotide gene sequences for the pea aphid using salmon version 1.1. For the 485 \npilot study, STAR (2.4.1b) [38] was used to map the reads to the pea aphid genome and HTseq counts 486 \nwas used to quantify the gene expression using AphidBase_OGS2.1b gene annotations.  487 \nClone-specific de novo  RNA-Seq assemblies (from both the heads and whole -body studies) were 488 \nindividually and collectively generated using Trinity version 2.9.1. All the data were pooled into one 489 \nfor the “collective” assembly, which was used for transcript differential expression analysis. The 490 \nindividual assemblies were used for gene prediction at a later stage. All RNA -Seq assemblies were 491 \nquality filtered using Transrate to reduce the probability of mis -assembled transcripts. Predicted 492 \ncoding sequences were generated using TransDecoder (with PFAM and BLAST guides). Diamond was 493 \nused to search against GenbankNR database. Differential expression analysis was performed using 494 \nEdgeR. Heatmaps and expression profile clustering w ere generated using the ptr script from within 495 \nthe Trinity package.  496 \nDuring early analysis, following visual assessment of RNA -seq read mapping and initial differential 497 \nexpression results, we found that the original pea aphid gene predictions (AphidBase_OGS2.1b) and 498 \nthe gene predictions fro m [39] did not fully match those generated by the de novo transcriptome 499 \nassemblies. Therefore, gene annotation was re -predicted on the published pea aphid genome 500 \n(OGS2.1b) to improve the accuracy of the gene models. Funannotate, in Other Eukaryotic mode, was 501 \nused to predict the genes using the de novo RNA-Seq assembly generated above, with RNA-Seq data 502 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n16 \n \nmapped using STAR (see above). A total of 29,930 genes were assigned codes in the format 503 \nACPISUM_0xxxxx, with the annotations provided at doi.org/10.5281/zenodo.11103500 [40]. 504 \nTo assign the various gene call s from the original genome assembly, bedtools intercept was used to 505 \nidentify genes with overlapping coordinates. If the genes overlapped, then they were considered the 506 \nsame gene. A simple BLAST appro ach could not be used here due to the duplicated nature of aphid 507 \nassemblies. A combination of reciprocal best BLAST hit, Orthofinder and MCL clustering were used to 508 \nassign genes between the clones as orthologues.  509 \nSaliva Collection 510 \nFor proteomics samples, N116 and PS01 were maintained separately on Vicia faba c.v. The Sutton, 511 \ngrown in standard potting compost and kept at 20 oC and a photoperiod of 16 -h light/8 -h dark. 512 \nApproximately 3,000 mixed aged aphids were positioned on 30 perspex rings (radius 4.5 cm, height 5 513 \ncm), each containing 4.5 ml of a chemically-defined diet, formulation A from [41], held between two 514 \nstretched sheets of ParafilmTM. The aphids were reared on the diets at 20°C with 18h light and 6h dark 515 \nfor 24 h after which  the diets were pooled and collected and stored at -80°C until required. Four 516 \nindependent replicates were produced by pooling the collected diet from two daily collections 517 \n(approximately 150 ml). Pooled diets were concentrated using a Vivacell 250 Pressure Concentrator 518 \n(Sartorius Mechatronics, UK) using a 5000 Da molecular weight cut -off (MWCO) polyethersulfone 519 \n(PES) membrane. When the final volume had reached 5 ml it was removed and 1 ml of filtered 520 \nsterilised PBS (phosphate-buffered saline) supplemented with Roche cOmplete TM protease inhibitor 521 \ncocktail (PIC) was added. The resulting mixture was further concentrated to approximately 250 μl 522 \nusing a Vivaspin 6 centrifuge concentrator (Sartorius Mechatronics, UK) with a 5000 Da MWCO PES 523 \nmembrane, purified using a 2D Clean -up Kit (GE HealthCare) following the manufacturer’s 524 \ninstructions. The resulting protein pellet was suspended in 25 μl 6 M urea, 2 M thiourea, 0.1 M Tris-525 \nHCl, pH 8.0 and re-quantified using the Qubit Fluorometer. Four independent biological replicates per 526 \ngenotype were subjected to mass spectrometry. 527 \nSalivary glands 528 \nThe salivary glands from 14-16 day old adult aphids of N116 and PS01 were dissected in ice-cold PBS 529 \nand transferred to 60 µl PBS supplemented with PIC. Forty pairs of salivary glands were pooled per 530 \nreplicate and homogenized with a motorised, disposable pestle. Sixty microliters of 12  M urea, 4 M 531 \nthiourea, and PIC was added and the samples were homogenised further and centrifuged at 9,000 × g 532 \nfor 5 min to pellet cellular debris. The supernatant was removed and quantified, and 100 µg of protein 533 \nwas purified using a 2D Clean -up Kit (GE HealthCare) following the manufacturer’s instructions with 534 \nthe exception that 400 μl of precipitant and co-precipitant were used in the first step . The resulting 535 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n17 \n \nprotein pellet was re-suspended in 30 μl 6 M urea, 2  M thiourea, 0.1 M Tris -HCl, pH 8.0 and re -536 \nquantified using the Qubit Fluorometer. Four biological replicate s per genotype were subjected to 537 \nmass spectrometry. 538 \nProtein sample digestion for mass spectrometry 539 \nThe digestion protocol was the same for both saliva and salivary gland samples and involved the 540 \naddition of 50 μl ammonium bicarbonate, reduction with 0.5 M dithiothreitol at 56°C for 20 min and 541 \nalkylation with 0.55 M iodoacetamide at room temperature for 15 min, in the dark. One μl of a 1% 542 \nw/v solution of ProteaseMax Surfactant Trypsin Enhancer (Promega) and 1 μg of Sequence Grade 543 \nTrypsin (Promega) were added , then samples were  incubated at 37°C for 18 h. Digestion was 544 \nterminated by adding 1 μl of 100% trichloroacetic acid (Sigma Aldrich) and incubati ng at room 545 \ntemperature for 5 min. Samples were centrifuged for 10 min at 13,000 x g and the supernatant was 546 \nremoved to new microcentrifuge tubes. 547 \nMass spectrometry and proteomic data analysis 548 \nOne μg of digested peptide was loaded onto a Dionex Ultimate 3000 (RSLCnano) chromatography 549 \nsystem connected to a QExactive (ThermoFisher Scientific) high -resolution accurate mass 550 \nspectrometer. Peptides were separated by an increasing acetonitrile gradient on a Biobasic C18 551 \nPicofritTM column (100 mm length, 75 µm ID), using 120 and 50 min reverse phase gradient s for 552 \nsalivary glands and saliva, respectively, at a flow rate of 250 nl min-1. All data were acquired with the 553 \nmass spectrometer operating in automatic data dependent switching mode. A high -resolution MS 554 \nscan (300 -2000 Da) was performed using the Orbitrap to select the 15 most intense ions prior to 555 \nMS/MS.  556 \nProtein identification and normalisation was conducted using the Andromeda search engine in 557 \nMaxQuant (version 1.6.17.0; http://maxquant.org/) to correlate the data against the predicted 558 \nprotein set generated in this study (ACPISUM_Proteins; 30891 entries) using default search 559 \nparameters for Orbitrap data.  False Discovery Rates were set to 1% for both peptides and proteins 560 \nand the FDR was estimated following searches against a target -decoy database. Two searches were 561 \nconducted for both N116 and PS01 saliva and salivary glands. The first involved a combined search of 562 \nthe raw files for each genotype separately to generate comprehensive proteomes for the saliva or 563 \nsalivary gland (hereafter All Identified Proteins). The second involved a quantitative search of the raw 564 \nfiles for all biological replicates (n=4) for the saliva or salivary glands. Quantitative and statistical 565 \nanalyses were conducted in Perseus (Version 1.6.1.1  http://maxquant.org/) using the n ormalized 566 \nlabel-free quantitation ( LFQ) intensity values from each sample . The data were filtered to remove 567 \ncontaminants, and peptides identified by site. LFQ intensity values were log2 transformed and samples 568 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n18 \n \nwere allocated to  their corresponding gr oups. A data imputation step was conducted to replace 569 \nmissing values with values that simulate signals of low abundant proteins chosen randomly from a 570 \ndistribution specified by a downshift of 2.1 times the mean standard deviation (SD) of all measured 571 \nvalues and a width of 0.1 times this SD. Normalized intensity values were used for principal 572 \ncomponents analysis. A two-sample t-test was performed using a cut-off value of p ≤ 0.05 to identify 573 \nstatistically significant differentially abundant (SSDA) proteins. Volcano plots were produced by 574 \nplotting –Log p-values on the y-axis and Log2 fold-change values on the x-axis to visualize differences 575 \nin protein abundance between the two genotypes.  576 \nGene annotations and Gene Ontology analysis  577 \nSecretion signal properties were predicted using SignalP4.1 [42]. Non-annotated genes were defined 578 \nas those with the following descriptors : hypothetical protein, uncharacterized protein , NA or 579 \nACYPIxxxxx without any other assigned function . GO enrichment analyses w ere performed using 580 \nGOseq [43]. 581 \nData availability 582 \nGenome annotations: zenodo.org/records/11103500 [40] 583 \nRNA-Seq: Pea aphid clones N116 and PS01 reared on Medicago truncatula A17 and DZA 315.16, 584 \ndissected heads: BioProject PRJNA757589, ncbi.nlm.nih.gov/bioproject/PRJNA757589/  585 \nRNA-Seq: Pea aphid clones N116 , PS01 and bulk F1 hybrid progeny reared on Medicago truncatula 586 \nA17, whole body samples: BioProject PRJNA757896, ncbi.nlm.nih.gov/bioproject/PRJNA757896  587 \nScripts: github.com/peterthorpe5/Pea_aphid_on_medicago_DZA_A17 588 \nProteomics: mass spectrometry data have been deposited to the ProteomeXchange Consortium via 589 \nthe PRIDE partner repository [44], dataset identifiers PXD053355 and PXD053620.  590 \n  591 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n19 \n \nReferences 592 \n1. Rossi M, Goggin FL, Milligan SB, Kaloshian I, Ullman DE, Williamson VM. The nematode resistance 593 \ngene Mi of tomato confers resistance against the potato aphid. Proc Natl Acad Sci. 1998;95:9750–4. 594 \n2. Dogimont C, Chovelon V, Pauquet J, Boualem A, Bendahmane A. The  V at  locus encodes for a CC ‐ 595 \nNBS ‐ LRR protein that confers resistance to  A phis gossypii  infestation and  A . gossypii  ‐mediated 596 \nvirus resistance. Plant J. 2014;80:993–1004. 597 \n3. Du B, Zhang W, Liu B, Hu J, Wei Z, Shi Z, et al. Identification and characterization of Bph14 , a gene 598 \nconferring resistance to brown planthopper in rice. Proc Natl Acad Sci. 2009;106:22163–8. 599 \n4. Rodriguez PA, Bos JIB. Toward Understanding the Role of Aphid Effectors in Plant Infestation. Mol 600 \nPlant-Microbe Interactions®. 2013;26:25–30. 601 \n5. Elzinga DA, De Vos M, Jander G. Suppression of Plant Defenses by a Myzus persicae (Green Peach 602 \nAphid) Salivary Effector Protein. Mol Plant-Microbe Interactions®. 2014;27:747–56. 603 \n6. Wang W, Dai H, Zhang Y, Chandrasekar R, Luo L, Hiromasa Y, et al. Armet is an effector protein 604 \nmediating aphid‐plant interactions. FASEB J. 2015;29:2032–45. 605 \n7. Mugford ST, Barclay E, Drurey C, Findlay KC, Hogenhout SA. An Immuno-Suppressive Aphid Saliva 606 \nProtein Is Delivered into the Cytosol of Plant Mesophyll Cells During Feeding. Mol Plant-Microbe 607 \nInteractions®. 2016;29:854–61. 608 \n8. Guo J, Wang H, Guan W, Guo Q, Wang J, Yang J, et al. A tripartite rheostat controls self-regulated 609 \nhost plant resistance to insects. Nature. 2023;618:799–807. 610 \n9. The International Aphid Genomics Consortium. Genome Sequence of the Pea Aphid 611 \nAcyrthosiphon pisum. PLoS Biol. 2010;8:e1000313. 612 \n10. Li Y, Park H, Smith TE, Moran NA. Gene Family Evolution in the Pea Aphid Based on 613 \nChromosome-Level Genome Assembly. Mol Biol Evol. 2019;36:2143–56. 614 \n11. Shigenobu S, Yorimoto S. Aphid hologenomics: current status and future challenges. Curr Opin 615 \nInsect Sci. 2022;50:100882. 616 \n12. Boulain H, Legeai F, Guy E, Morlière S, Douglas NE, Oh J, et al. Fast Evolution and Lineage-Specific 617 \nGene Family Expansions of Aphid Salivary Effectors Driven by Interactions with Host-Plants. Genome 618 \nBiol Evol. 2018;10:1554–72. 619 \n13. Boulain H, Legeai F, Jaquiéry J, Guy E, Morlière S, Simon J-C, et al. Differential Expression of 620 \nCandidate Salivary Effector Genes in Pea Aphid Biotypes With Distinct Host Plant Specificity. Front 621 \nPlant Sci. 2019;10:1301. 622 \n14. Wang D, Yang Q, Hu X, Liu B, Wang Y. A Method for Identification of Biotype-Specific Salivary 623 \nEffector Candidates of Aphid. Insects. 2023;14:760. 624 \n15. Bos JIB, Prince D, Pitino M, Maffei ME, Win J, Hogenhout SA. A Functional Genomics Approach 625 \nIdentifies Candidate Effectors from the Aphid Species Myzus persicae (Green Peach Aphid). PLOS 626 \nGenet. 2010;6:e1001216. 627 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n20 \n \n16. Thorpe P, Cock PJA, Bos J. Comparative transcriptomics and proteomics of three different aphid 628 \nspecies identifies core and diverse effector sets. BMC Genomics. 2016;17:172. 629 \n17. Nicolis VF, Burger NFV, Botha A-M. Whole-body transcriptome mining for candidate effectors 630 \nfrom Diuraphis noxia. BMC Genomics. 2022;23:493. 631 \n18. Rotenberg D, Baumann AA, Ben-Mahmoud S, Christiaens O, Dermauw W, Ioannidis P, et al. 632 \nGenome-enabled insights into the biology of thrips as crop pests. BMC Biol. 2020;18:142. 633 \n19. Grbić M, Van Leeuwen T, Clark RM, Rombauts S, Rouzé P, Grbić V, et al. The genome of 634 \nTetranychus urticae reveals herbivorous pest adaptations. Nature. 2011;479:487–92. 635 \n20. Huo S-M, Yan Z-C, Zhang F, Chen L, Sun J-T, Hoffmann AA, et al. Comparative genome and 636 \ntranscriptome analyses reveal innate differences in response to host plants by two color forms of the 637 \ntwo-spotted spider mite Tetranychus urticae. BMC Genomics. 2021;22:569. 638 \n21. Mathers TC, Chen Y, Kaithakottil G, Legeai F, Mugford ST, Baa-Puyoulet P, et al. Rapid 639 \ntranscriptional plasticity of duplicated gene clusters enables a clonally reproducing aphid to colonise 640 \ndiverse plant species. Genome Biol. 2017;18:27. 641 \n22. Peccoud J, Ollivier A, Plantegenest M, Simon J-C. A continuum of genetic divergence from 642 \nsympatric host races to species in the pea aphid complex. Proc Natl Acad Sci. 2009;106:7495–500. 643 \n23. Caillaud MC, Mondor‐Genson G, Levine‐Wilkinson S, Mieuzet L, Frantz A, Simon JC, et al. 644 \nMicrosatellite DNA markers for the pea aphid Acyrthosiphon pisum. Mol Ecol Notes. 2004;4:446–8. 645 \n24. Eyres I, Jaquiéry J, Sugio A, Duvaux L, Gharbi K, Zhou J, et al. Differential gene expression 646 \naccording to race and host plant in the pea aphid. Mol Ecol. 2016;25:4197–215. 647 \n25. Eyres I, Duvaux L, Gharbi K, Tucker R, Hopkins D, Simon J, et al. Targeted re‐sequencing confirms 648 \nthe importance of chemosensory genes in aphid host race differentiation. Mol Ecol. 2017;26:43–58. 649 \n26. Kanvil S, Powell G, Turnbull C. Pea aphid biotype performance on diverse Medicago host 650 \ngenotypes indicates highly specific virulence and resistance functions. Bull Entomol Res. 651 \n2014;104:689–701. 652 \n27. Kanvil S, Collins CM, Powell G, Turnbull CGN. Cryptic Virulence and Avirulence Alleles Revealed 653 \nby Controlled Sexual Recombination in Pea Aphids. Genetics. 2015;199:581–93. 654 \n28. Stewart SA, Hodge S, Ismail N, Mansfield JW, Feys BJ, Prospéri J-M, et al. The RAP1 Gene Confers 655 \nEffective, Race-Specific Resistance to the Pea Aphid in Medicago truncatula Independent of the 656 \nHypersensitive Reaction. Mol Plant-Microbe Interactions®. 2009;22:1645–55. 657 \n29. Julier B, Huguet T, Chardon F, Ayadi R, Pierre J-B, Prosperi J-M, et al. Identification of quantitative 658 \ntrait loci influencing aerial morphogenesis in the model legume Medicago truncatula. Theor Appl 659 \nGenet. 2007;114:1391–406. 660 \n30. Peccoud J, Mahéo F, de la Huerta M, Laurence C, Simon J-C. Genetic characterisation of new 661 \nhost-specialised biotypes and novel associations with bacterial symbionts in the pea aphid complex. 662 \nInsect Conserv Divers. 2015;8:484–92. 663 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n21 \n \n31. Carolan JC, Caragea D, Reardon KT, Mutti NS, Dittmer N, Pappan K, et al. Predicted Effector 664 \nMolecules in the Salivary Secretome of the Pea Aphid ( Acyrthosiphon pisum ): A Dual 665 \nTranscriptomic/Proteomic Approach. J Proteome Res. 2011;10:1505–18. 666 \n32. Carolan JC, Fitzroy CIJ, Ashton PD, Douglas AE, Wilkinson TL. The secreted salivary proteome of 667 \nthe pea aphid  Acyrthosiphon pisum  characterised by mass spectrometry. PROTEOMICS. 668 \n2009;9:2457–67. 669 \n33. Pitino M, Hogenhout SA. Aphid Protein Effectors Promote Aphid Colonization in a Plant Species-670 \nSpecific Manner. Mol Plant-Microbe Interactions®. 2013;26:130–9. 671 \n34. Kaloshian I, Walling LL. Hemipterans as Plant Pathogens. Annu Rev Phytopathol. 2005;43:491–672 \n521. 673 \n35. Furch ACU, Van Bel AJE, Will T. Aphid salivary proteases are capable of degrading sieve-tube 674 \nproteins. J Exp Bot. 2015;66:533–9. 675 \n36. Cantón PE, Bonning BC. Extraoral digestion: outsourcing the role of the hemipteran midgut. Curr 676 \nOpin Insect Sci. 2020;41:86–91. 677 \n37. Ferrari J, Via S, Godfray HCJ. POPULATION DIFFERENTIATION AND GENETIC VARIATION IN 678 \nPERFORMANCE ON EIGHT HOSTS IN THE PEA APHID COMPLEX. Evolution. 2008;62:2508–24. 679 \n38. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-680 \nseq aligner. Bioinformatics. 2013;29:15–21. 681 \n39. Thorpe P, Escudero-Martinez CM, Cock PJA, Eves-van Den Akker S, Bos JIB. Shared 682 \nTranscriptional Control and Disparate Gain and Loss of Aphid Parasitism Genes. Genome Biol Evol. 683 \n2018;10:2716–33. 684 \n40. Thorpe P. Genome annotations for: Multi-omics approaches define novel aphid effector 685 \ncandidates associated with virulence and avirulence phenotypes. 2024. 686 \n41. Prosser WA, Douglas AE. The aposymbiotic aphid: An analysis of chlortetracycline-treated pea 687 \naphid, Acyrthosiphon pisum. J Insect Physiol. 1991;37:713–9. 688 \n42. Petersen TN, Brunak S, von Heijne G, Nielsen H. SignalP 4.0: discriminating signal peptides from 689 \ntransmembrane regions. Nat Methods. 2011;8:785–6. 690 \n43. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting 691 \nfor selection bias. Genome Biol. 2010;11:R14. 692 \n44. Perez-Riverol Y, Bai J, Bandla C, García-Seisdedos D, Hewapathirana S, Kamatchinathan S, et al. 693 \nThe PRIDE database resources in 2022: a hub for mass spectrometry-based proteomics evidences. 694 \nNucleic Acids Res. 2022;50:D543–52. 695 \n  696 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n22 \n \nTables 697 \nTable 1 Genes and proteins overlapping in multiple experiments. All genes shown that are 698 \nrepresented in at least three datasets, plus all genes intersected between F1 transcriptome and at 699 \nleast one other dataset. Saliva and salivary gland data are proteins, head and body data are 700 \ntranscripts. A. Proteins and upregulated genes in virulent aphids (N116, VIR F1 pool); B. Proteins and 701 \nupregulated genes in avirulent aphids (PS01, AVR F1 pool). Y = protein present and/or RNA 702 \ndifferentially expressed. Full gene and protein lists are in Supplementary Material 3 and 4. 703 \nGene Annotation Saliva Salivary \ngland \nParent \nhead \nParent \nbody \nF1 \nbody \nA) N116 & VIR F1       \nACPISUM_000319 ACYPI007553 Y Y  Y  \nACPISUM_006458 aldo-keto reductase family 1 member B10-like  Y Y Y  \nACPISUM_025240 aminopeptidase N Y Y Y Y  \nACPISUM_005699 aminopeptidase N Y Y Y Y  \nACPISUM_025168 aminopeptidase N Y Y Y Y  \nACPISUM_009258 aminopeptidase N Y  Y Y  \nACPISUM_024778 aminopeptidase N Y  Y Y  \nACPISUM_026844 aminopeptidase N Y  Y Y  \nACPISUM_025015 aminopeptidase N Y Y Y Y  \nACPISUM_023906 Apoptosis inducing protein  Y Y Y  \nACPISUM_020864 F-actin-capping protein subunit alpha  Y Y Y  \nACPISUM_023535 glutamate-gated chloride channel-like   Y  Y \nACPISUM_019971 glutathione hydrolase 1 proenzyme-like  Y Y Y  \nACPISUM_010531 hypothetical protein X975_16721    Y Y \nACPISUM_013751 LYR motif-containing protein 4  Y Y Y  \nACPISUM_013796 myrosinase 1-like    Y Y \nACPISUM_006164 ---NA---  Y Y Y  \nACPISUM_023321 papain inhibitor-like  Y Y Y  \nACPISUM_009624 proline-rich extensin-like protein EPR1  Y Y Y  \nACPISUM_028519 single-stranded DNA-binding replication protein A  Y Y Y  \nACPISUM_025560 ubiquinone biosynthesis monooxygenase COQ6, \nmitochondrial \n Y Y Y  \nACPISUM_008675 uncharacterized protein LOC100162547   Y Y Y  \nACPISUM_007320 uncharacterized protein LOC100167449  Y Y Y Y  \nACPISUM_001031 uncharacterized protein LOC100571631   Y Y Y  \nACPISUM_016519 uncharacterized protein LOC100573156  Y Y Y  \nACPISUM_010687 uncharacterized protein LOC103309122  Y Y Y  \nACPISUM_017388 uncharacterized protein LOC103309964 Y Y Y Y  \nACPISUM_009099 uncharacterized protein LOC112598674  Y Y Y  \nACPISUM_027918 vacuolar protein sorting-associated protein 29  Y Y Y  \nB) PS01 & AVR F1       \nACPISUM_000957 AGAP002382-PA-like protein   Y Y Y \nACPISUM_015173 AGAP011571-PA-like protein   Y  Y \nACPISUM_002223 aminopeptidase N  Y Y Y  \nACPISUM_003737 aminopeptidase N Y Y Y Y  \nACPISUM_023448 aminopeptidase N Y Y Y Y  \nACPISUM_028967 aminopeptidase N Y Y Y Y  \nACPISUM_021545 aminopeptidase N Y Y Y   \nACPISUM_009259 aminopeptidase N  Y Y Y  \nACPISUM_009580 anoctamin-1-like   Y Y Y \nACPISUM_012705 CD63 antigen   Y Y Y \nACPISUM_006933 cuticular protein   Y Y Y \nACPISUM_019160 glutathione S-transferase 1-1-like Y Y Y Y  \nACPISUM_019168 glutathione S-transferase 1-1-like Y Y Y Y  \nACPISUM_001883 glutathione S-transferase D7-like Y Y Y   \nACPISUM_016389 histone acetyltransferase KAT6B isoform X1    Y Y \nACPISUM_009097 multidrug resistance-associated protein 1  Y Y Y  \nACPISUM_011553 ---NA---  Y Y Y  \nACPISUM_011754 ---NA--- Y Y Y Y  \nACPISUM_004702 ---NA---    Y Y \nACPISUM_021569 ---NA---    Y Y \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n23 \n \nACPISUM_025236 ---NA---    Y Y \nACPISUM_014327 ---NA---   Y Y Y \nACPISUM_016390 ---NA---   Y Y Y \nACPISUM_017200 ---NA---   Y Y Y \nACPISUM_027631 ---NA---   Y Y Y \nACPISUM_028853 ---NA---   Y Y Y \nACPISUM_019381 neural cell adhesion molecule L1 isoform X1   Y Y Y \nACPISUM_020816 peroxidase-like  Y Y Y Y  \nACPISUM_019857 peroxidase-like Y Y Y   \nACPISUM_019870 peroxidasin homolog  Y Y Y  \nACPISUM_000958 phospholipase DDHD2-like    Y  Y \nACPISUM_006758 piggyBac transposable element-derived protein 4-like   Y Y Y \nACPISUM_022113 piwi-like protein Siwi    Y Y \nACPISUM_010778 predicted protein   Y Y Y \nACPISUM_019013 protein ABHD18   Y Y Y \nACPISUM_021997 regucalcin-like  Y Y  Y Y \nACPISUM_021999 regucalcin-like Y Y Y   \nACPISUM_001383 replication protein A 70 kDa DNA-binding subunit-like  Y Y Y  \nACPISUM_015166 TBC1 domain family member 19   Y Y Y \nACPISUM_014232 tubulin glycylase 3A-like    Y Y \nACPISUM_008377 uncharacterized family 31 glucosidase KIAA1161-like   Y Y Y \nACPISUM_008379 uncharacterized family 31 glucosidase KIAA1161-like   Y Y Y \nACPISUM_008380 uncharacterized family 31 glucosidase KIAA1161-like   Y Y Y \nACPISUM_012348 uncharacterized protein LOC100158692 Y Y Y   \nACPISUM_018433 uncharacterized protein LOC100158721    Y Y \nACPISUM_007487 uncharacterized protein LOC100160601 Y Y Y   \nACPISUM_016065 uncharacterized protein LOC100161530  Y Y Y  \nACPISUM_007076 uncharacterized protein LOC100163035   Y Y Y \nACPISUM_016064 uncharacterized protein LOC100570074  Y Y Y  \nACPISUM_029311 uncharacterized protein LOC100570454 Y Y Y Y  \nACPISUM_008664 uncharacterized protein LOC100570454 Y Y Y   \nACPISUM_007394 uncharacterized protein LOC100572241 Y Y Y Y  \nACPISUM_021703 uncharacterized protein LOC100575642 Y  Y Y  \nACPISUM_029930 uncharacterized protein LOC100575698 Y Y Y Y Y \nACPISUM_006906 uncharacterized protein LOC100575848   Y Y Y \nACPISUM_003989 uncharacterized protein LOC103307823  Y Y Y  \nACPISUM_024374 uncharacterized protein LOC107882950    Y Y \nACPISUM_015285 uncharacterized protein LOC107883982   Y Y Y \nACPISUM_000491 uncharacterized protein LOC111028731  Y Y Y  \nACPISUM_027814 uncharacterized SDCCAG3 family protein-like   Y Y Y \n 704 \n  705 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n24 \n \nTable 2 Significantly enriched GO terms within differentially expressed transcript and 706 \nprotein data. Terms enriched at FDR<0.05, after manual curation to remove redundancies, 707 \nretaining the terms with lowest FDR. Full lists of enriched terms are in Supplementary 708 \nMaterial 5. 709 \nGO Category Term \nOntology \ngroup \nNo. genes \nin DE set P value FDR \nWhole body RNA      \nN116 up      \n0004177 aminopeptidase activity MF 8 2.72E-06 0.0270 \nPS01 up      \n0004177 aminopeptidase activity MF 9 1.43E-07 0.0014 \n0017177 glucosidase II complex CC 4 2.09E-06 0.0104 \nAVR F1 up      \n0017177 glucosidase II complex CC 3 2.42E-06 0.0240 \nVIR F1 up No terms     \nHeads RNA      \nN116 up on A17      \n0045271 respiratory chain complex I CC 11 1.55E-07 9.63E-05 \n0005743 mitochondrial inner membrane CC 22 2.01E-07 0.0001 \n0004177 aminopeptidase activity MF 10 2.13E-06 0.0009 \n0016491 oxidoreductase activity MF 31 3.24E-06 0.0012 \n0005875 microtubule associated complex CC 25 1.72E-05 0.0052 \n0019395 fatty acid oxidation BP 5 2.44E-05 0.0064 \n0042826 histone deacetylase binding MF 4 0.00011 0.0275 \n0045239 tricarboxylic acid cycle enzyme complex CC 3 0.00014 0.0327 \n0004448 isocitrate dehydrogenase activity MF 3 0.00015 0.0338 \nN116 up on DZA      \n0006635 fatty acid beta-oxidation BP 5 1.05E-06 0.0082 \n0004177 aminopeptidase activity MF 9 2.66E-06 0.0082 \n0004449 isocitrate dehydrogenase (NAD+) activity MF 3 1.40E-05 0.0198 \n0006099 tricarboxylic acid cycle BP 6 3.54E-05 0.0389 \nPS01 up on A17      \n0004177 aminopeptidase activity MF 10 7.77E-07 0.00771 \nPS01 up on DZA      \n0004177 aminopeptidase activity MF 11 5.59E-09 5.53E-05 \nSalivary gland proteins      \nN116 up      \n0003983 UTP:glucose-1-phosphate \nuridylyltransferase activity \nMF 3 4.97E-06 0.0213 \nPS01 up No terms     \n 710 \n  711 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n25 \n \nTable 3. Comparison of expression patterns of exopeptidases detected in saliva and salivary 712 \nglands. All detected proteins are listed, along with whether they were differentially expressed, and 713 \nwhether the patterns were also reflected in the transcriptomes. Sal = saliva; SG = salivary gland; Y = 714 \nprotein present. 715 \n 716 \n  717 \nACPISUM v3 gene Nearest ACYPI gene(s) Sal Sal SG SG Sal SG heads whole\nN116 PS01 N116 PS01\nAminopeptidase N\nACPISUM_005699-T1 ACYPI080623 ACYPI070600 ACYPI005810 Y Y Y Y\nACPISUM_025168-T1 ACYPI068031 Y Y\nACPISUM_027632-T1 ACYPI073645 Y Y\nACPISUM_009259-T1 ACYPI007868 Y Y Y\nACPISUM_009258-T1 ACYPI007868 Y Y\nACPISUM_024778-T1 ACYPI072916 Y\nACPISUM_025015-T1 ACYPI061522 ACYPI21510 Y Y\nACPISUM_023174-T1 ACYPI006366 Y Y\nACPISUM_025240-T1 ACYPI070333 Y Y Y\nACPISUM_000115-T1 ACYPI000001 Y\nACPISUM_003737-T1 ACYPI067691 Y Y Y\nACPISUM_021545-T1 ACYPI085147 ACYPI002583 Y Y Y Y\nACPISUM_029674-T1 ACYPI072988 Y Y Y Y\nACPISUM_023448-T1 ACYPI010198 Y Y Y Y\nACPISUM_028967-T1 ACYPI083965 Y Y Y Y\nACPISUM_000246-T1 ACYPI086097 ACYPI43770 ACYPI068046 Y Y\nACPISUM_010796-T1 ACYPI071232 ACYPI33244 Y Y Y Y\nACPISUM_012062-T1 ACYPI22813 Y Y Y Y\nACPISUM_018507-T1 ACYPI21711 ACYPI084528 ACYPI003165 Y Y Y Y\nACPISUM_006298-T1 ACYPI41708 ACYPI22605 Y Y Y Y\nACPISUM_019635-T1 ACYPI060722 Y\nACPISUM_019937-T1 ACYPI49161 Y Y Y Y\nACPISUM_026119-T1 ACYPI083984 Y Y Y Y\nACPISUM_017858-T1 ACYPI54528 ACYPI001911 Y Y\nACPISUM_002219-T1 ACYPI44040 Y Y Y Y\nACPISUM_014203-T1 ACYPI067721 Y Y\nACPISUM_018506-T1 ACYPI21557 Y Y\nACPISUM_019609-T1 ACYPI001203 Y Y\nACPISUM_019610-T1 ACYPI071951 Y Y\nAminopeptidase N total proteins detected 24 20 21 17\nAngiotensin converting enzyme\nACPISUM_008374-T1 ACYPI000733 Y Y Y Y\nACPISUM_024301-T1 ACYPI084554 Y Y Y Y\nACPISUM_024303-T1 ACYPI071320 Y Y Y Y\nACPISUM_020790-T1 ACYPI008911 Y Y\nKey\nnot detected\ndetected, not DE\nup in N116\nup in PS01\nPresence/absence Differential expression\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n26 \n \nFigure legends 718 \nFigure 1. Summary of transcriptome and proteome analysis pipeline. Resistant and susceptible host 719 \nplants carried or lacked the RAP1 aphid resistance QTL, respectively. For all experiments, virulent N116 720 \nand avirulent PS01 aphids were compared. In addition, BSA-RNA-Seq was done on whole body pooled 721 \nsamples of F1 virulent and avirulent aphids. 722 \nFigure 2. Virulence phenotypes of parental aphid clones and selections from F1 populations used 723 \nfor BSA-RNA-Seq. Tested on two M. truncatula genotypes carrying the RAP1 locus: Jemalong A17 and 724 \na resistant near isogenic line (RNIL) derived from a cross between A17 and DZA315.16. The parental 725 \ngenotypes and selections from the F1 populations shown here were all used for the BSA -RNA-Seq 726 \nexperiment. Data are expressed as virulence index, assessed 10 d after infestation. Phenotypes of F1 727 \nclones were classified using the following virulence index cut -offs: A17 VIR >4, AVR <2; RNIL VIR >4, 728 \nAVR <4.5. Orange circles are NP (N116 female x PS01 male); blue triangles are PN (PS01 female x N116 729 \nmale); red is N116, and green is PS01, with each of three parental data points from a separate batch 730 \nof F1 tests. The full population phenotype data are provided in Supplementary Material 1. 731 \nFigure 3. Transcriptome analysis of aphid heads. Samples were dissected heads from PS01 and N116 732 \ngenotypes infested on Medicago truncatula A17 or DZA315.16 for 24 h, with n=3 biological replicates. 733 \nAphid genotype PS01 is avirulent on M. truncatula  A17, and all other combinations represent 734 \ncompatible interactions. A. Clustering of transcriptional responses of pea aphid, showing samples 735 \nclustered more strongly based on aphid genotype than on host interaction; B. Principal components 736 \nanalysis. The top two principal components explain >68% of the variation among transcriptional 737 \nresponses. Samples group largely by aphid genotype rather than host interaction; C. Numbers of genes 738 \ndifferentially expressed between the different aphid genotypes on different M. truncatula genotypes. 739 \nOf the 935 DE genes between PS01 and N116 on A17, 483 were up in N116 and 452 were up in PS01. 740 \nOf the 758 DE genes on DZA hosts, 395 were up in N116 and 363 were up in PS01. Accompanying gene 741 \nlists and annotations are provided in Supplementary Material 3. 742 \nFigure 4. Transcriptome analysis of whole aphids. Aphids were infested on Medicago truncatula A17 743 \nfor 24 h, with n=5 biological replicates. A . Clustering of transcriptional responses of pea aphid PS01, 744 \nN116, bulked F1 VIR and AVR progeny. Responses within biological replicates are more strongly 745 \ncorrelated than responses among different aphid genotypes; B. Principal components analysis. The 746 \ntop 2 principal components explain >45% of the variation among transcriptional responses of PS01, 747 \nN116, AVR and VIR F1 progeny replicates, and separate the responses of the different aphid genotypes 748 \nand F1 pools. C. Numbers of genes differentially expressed between the different aphid genotypes 749 \nand pools. Accompanying gene lists and annotations are provided in Supplementary Material 3. 750 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n27 \n \nFigure 5. Differential gene expression in pea aphid genotypes N116, PS01 and bulked F1 pools of 751 \nvirulent and avirulent progeny. A. Numbers of genes up - versus down-regulated in comparisons of 752 \nparent genotypes N116 (virulent) and PS01 (avirulent), and their F1 progeny pools, all on A17 host 753 \nplants. Orange bars represent the numbers of genes up-regulated in genotype N116 or the VIR F1 pool 754 \ncompared with the genotype PS01 and the AVR pool, respectively. Blue bars represent the numbers 755 \nof genes up-regulated in genotype PS01 or the AVR F1 pool compared with the genotype N116 and 756 \nthe VIR pool, respectively; B. Overlaps in genes up-regulated in genotype N116 whole body and head 757 \ntissues compared to genotype PS01, and up-regulated in the VIR F1 pool compared to the AVR F1 pool; 758 \nC. Overlaps in genes up -regulated in genotype PS01 whole body and head tissues compared to 759 \ngenotype N116, and up-regulated in the AVR F1 pool compared to the VIR F1 pool; D. Overlaps in up-760 \nregulated genes among whole body transcriptomes of N116, PS01, VIR F1 pool and AVR F1 pool. 761 \nFigure 6. Comparative proteomic analysis of salivary glands and saliva for pea aphid genotypes N116 762 \nand PS01. Venn diagrams of the number of proteins shared and found exclusively in A) salivary glands 763 \nand B) saliva identified for both genotypes. Principal Components Analysis (PCA) of C) salivary glands 764 \nand D) saliva distinguishes both genotypes clearly. Volcano plots based on -log10 p values and log2 fold 765 \ndifferences highlighting the statistically significant differentially abundant (SSDA) proteins (p≤0.05) for 766 \nE) salivary glands and F) saliva. Annotations are shown for the top 12 proteins of increased and 767 \ndecreased abundances. 768 \nFigure 7. Selected differentially expressed genes from whole body transcriptomes . A,B 769 \nrepresentative genes upregulated both in virulent parent and in virulent F1 pool; C -F representative 770 \ngenes upregulated both in avirulent parent and in avirulent F1 pool. G,H representative genes with 771 \nopposite regulation between parent and F1 pairs. Each point represents an individual RNA-Seq library 772 \n(n=5). *** indicates FDR<0.001. 773 \n  774 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n28 \n \nAcknowledgements  775 \nFunding 776 \nWe thank the Biotechnology and Biological Sciences Research Council for funding to CT 777 \n(BB/N002830/1) and JB (BB/N002660/1). We thank Umer Rashid and Martin Selby for expert 778 \ntechnical assistance. 779 \nAuthor contributions 780 \nJB, CT, JC, PT, SA and RLC designed the experiments. PT, SA, RLC, ND, JCS, JI and SK conducted the 781 \nexperiments. PT, SA, RLC, JB, CT, JC, ND and JCS analysed the data. CT, JB, JC and PT wrote the paper. 782 \nAll authors approved the submitted manuscript. 783 \nConflicts 784 \nThe authors declare that they have no competing interests. 785 \nSupplementary Materials 786 \nSupplementary Material 1. F1 aphid phenotyping. 787 \nSupplementary Material 2. Read mapping summary. 788 \nSupplementary Material 3. RNA-Seq Head and whole body differentially expressed genes. 789 \nSupplementary Material 4. Salivary gland and saliva proteomics. 790 \nSupplementary Material 5. GO enrichment791 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n29 \n \n 792 \n 793 \nFigure 1. Summary of transcriptome and proteome analysis pipeline. Resistant and susceptible host plants carried or lacked the RAP1 aphid resistance QTL, 794 \nrespectively. For all experiments, virulent N116 and avirulent PS01 aphids were compared. In addition, BSA-RNA-Seq was done on whole body pooled samples 795 \nof F1 virulent and avirulent aphids. 796 \n  797 \nVirulent vs\navirulent aphids\nDifferential protein \nabundance\nProteomics\nHeads\nRNA-Seq\nSalivary \ngland SalivaWhole \nbody\nDifferential gene \nexpression \nResistant vs \nsusceptible host\nCandidate Effectors\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n30 \n \n 798 \n 799 \nFigure 2. Virulence phenotypes of parental aphid clones and selections from F1 populations used for BSA-RNA-Seq. Tested on two M. truncatula genotypes 800 \ncarrying the RAP1 locus: Jemalong A17 and a resistant near isogenic line (RNIL) derived from a cross between A17 and DZA315.16. The parental ge notypes 801 \nand selections from the F1 populations shown here were all used for the BSA-RNA-Seq experiment. Data are expressed as virulence index, assessed 10 d after 802 \ninfestation. Phenotypes of F1 clones were classified using the following virulence index cut -offs: A17 VIR >4, AVR <2; RNIL VIR >4, AVR <4.5. Orange circles 803 \nare NP (N116 female x PS01 male); blue triangles are PN (PS01 female x N116 male); red is N116, and green is PS01, with each of three parental data points 804 \nfrom a separate batch of F1 tests. The full population phenotype data are provided in Supplementary Material 1. 805 \n  806 \n0\n1\n2\n3\n4\n5\n6\n \n8\n9\n0 1 2 3 4 5 6  8 9\nVirulence index on rNIL\nVirulence index on A1 \n    \n    \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n31 \n \n 807 \nFigure 3. Transcriptome analysis of aphid heads. Samples were dissected heads from PS01 and N116 genotypes infested on Medicago truncatula A17 or 808 \nDZA315.16 for 24 h, with n=3 biological replicates. Aphid genotype PS01 is avirulent on M. truncatula A17, and all other combinations represent compatible 809 \ninteractions. A. Clustering of transcriptional responses of pea aphid, showing samples clustered more strongly based on aphid  genotype than on host 810 \ninteraction; B. Principal components analysis. The top two principal components explain >68% of the variation among transcrip tional responses. Samples 811 \ngroup largely by aphid genotype rather than host interaction; C. Numbers of genes differentially expressed between the different aphid genotypes on different 812 \nM. truncatula genotypes. Of the 935 DE genes between PS01 and N116 on A17, 483 were up in N116 and 452 were up in PS01. Of the 758 DE genes  on DZA 813 \nhosts, 395 were up in N116 and 363 were up in PS01. Accompanying gene lists and annotations are provided in Supplementary Material 3.   814 \nA B\nPS01/DZAPS01/A17N116/DZAN116/A17\n845935330N116/A17\n758782033N116/DZA\n290782935PS01/A17\n029758845PS01/DZA\nC\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n32 \n \n 815 \nFigure 4. Transcriptome analysis of whole aphids. Aphids were infested on Medicago truncatula A17 for 24 h, with n=5 biological replicates. A. Clustering of 816 \ntranscriptional responses of pea aphid PS01, N116, bulked F1 VIR and AVR progeny. Responses within biological replicates are more strongly correlated than 817 \nresponses among different aphid genotypes; B. Principal components analysis. The top 2 principal components explain >45% of t he variation among 818 \ntranscriptional responses of PS01, N116, AVR and VIR F1 progeny replicates, and separate the responses of the different aphid  genotypes and F1 pools. C. 819 \nNumbers of genes differentially expressed between the different aphid genotypes and pools. Accompanying gene lists and annota tions are provided in 820 \nSupplementary Material 3.  821 \n822 \nB\nAvr F1 \npoolPS01Vir F1 \npoolN116\n2695771480N116\n882410148Vir F1 pool\n2460241577PS01\n024688269Avr F1 pool\nC\nA\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n33 \n \n 823 \nFigure 5. Differential gene expression in pea aphid genotypes N116, PS01 and bulked F1 pools of virulent and avirulent progeny. A. Numbers of genes up- 824 \nversus down-regulated in comparisons of parent genotypes N116 (virulent) and PS01 (avirulent), and their F1 progeny pools, all on A17 hos t plants. Orange 825 \nbars represent the numbers of genes up -regulated in genotype N116 or the VIR F1 pool compared with the genotype PS01 and the AVR pool, respectively. 826 \nBlue bars represent the numbers of genes up-regulated in genotype PS01 or the AVR F1 pool compared with the genotype N116 and the VIR pool, respectively; 827 \nB. Overlaps in genes up-regulated in genotype N116 whole body and head tissues compared to genotype PS01, and up-regulated in the VIR F1 pool compared 828 \nto the AVR F1 pool; C. Overlaps in genes up -regulated in genotype PS01 whole body and head tissues compared to genotype N116, and up -regulated in the 829 \nAVR F1 pool compared to the VIR F1 pool; D. Overlaps in up-regulated genes among whole body transcriptomes of N116, PS01, VIR F1 pool and AVR F1 pool. 830 \n132\n165\n317\n21\n2\n1\nN116 up - heads\nN116 up - whole aphid\nVir F1 pool - up\nB\n260\n168\n79\n31\n9\n22\n2\nPS01 up - whole aphid\nPS01 up - heads\nAvr F1 pool - up\nC\n-500\n-400\n-300\n-200\n-100\n0\n100\n200\n300\n400\n500\nA\n299\n24\n452\n64\nnumber of DE genes\navirulentvirulent\n278\n483\nD N116 up - whole aphid\nPS01 up - whole aphid\nAvr F1 pool - up\nVir F1 pool -up\n294\n234\n31\n139\n30\n2\n3\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n34 \n \n 831 \nFigure 6. Comparative proteomic analysis of salivary glands and saliva for pea aphid genotypes N116 and PS01. Venn diagrams of the number of proteins shared and found exclusively in A) salivary glands and B) saliva identified for 832 \nboth genotypes. Principal Components Analysis ( PCA) of C) salivary glands and D) saliva distinguishes both genotypes clearly. Volcano plots based on -log10 p values and log 2 fold differences highlighting the statistically significant 833 \ndifferentially abundant (SSDA) proteins (p≤0.05) for E) salivary glands and F) saliva. Annotations are shown for the top 12 proteins of increased and decreased abundances.834 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint \n\n35 \n \n 835 \nFigure 7. Selected differentially expressed genes from whole body transcriptomes . A,B representative genes upregulated both in virulent parent and in 836 \nvirulent F1 pool; C-F representative genes upregulated both in avirulent parent and in avirulent F1 pool. G,H representative genes with opposite regulation 837 \nbetween parent and F1 pairs. Each point represents an individual RNA-Seq library (n=5). *** indicates FDR<0.001. 838 \n‐2\n‐1\n0\n1\n2\n3\n ero centred log2 \nexpression\nN116        PS01      F1 VIR   F1 AVR\nACPIS M 006933\ncu cular protein\n              \n‐3\n‐2\n‐1\n0\n1\n2\n3\n ero centred log2 \nexpression\nN116        PS01      F1 VIR   F1 AVR\nACPIS M 0198 0\nperoxidasin homolog\n              \n‐2\n‐1\n0\n1\n2\n ero centred log2 \nexpression\nN116        PS01      F1 VIR   F1 AVR\nACPIS M 005464\ntranscrip onal regulator ERG homolog\n              \n‐3\n‐2\n‐1\n0\n1\n2\n3\n ero centred log2 \nexpression\nN116        PS01      F1 VIR   F1 AVR\nACPIS M 013 96\nmyrosinase 1‐like\n              \n‐2\n‐1\n0\n1\n2\n ero centred log2 \nexpression\nN116        PS01      F1 VIR   F1 AVR\nACPIS M 0199 1\nglutathione hydrolase 1 proenzyme‐like \n              \n‐4\n‐3\n‐2\n‐1\n0\n1\n2\n3\n ero centred log2 \nexpression\nN116        PS01      F1 VIR   F1 AVR\nACPIS M 008380\nfamily 31 glucosidase KIAA1161‐like\n              \n‐2\n‐1\n0\n1\n2\n ero centred log2 \nexpression\nN116        PS01      F1 VIR   F1 AVR\nACPIS M Api15 8  \nuncharacterisedprotein\n              \n‐3\n‐2\n‐1\n0\n1\n2\n3\n ero centred log2 \nexpression\nN116        PS01      F1 VIR   F1 AVR\nACPIS M 02199 \nregucalcin‐like\n              \nA B\nC D\nE F\nG H\nACPISUM_029930\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}