Multi-omics approaches define novel aphid effector candidates associated with virulence and avirulence phenotypes

12 Background. Compatibility between plant parasites and their hosts is genetically determined by both 13 interacting organisms. For example, plants may carry resistance (R) genes or deploy chemical 14 defences. Aphid saliva contains many proteins that are secreted into host tissues. Subsets of these 15 proteins are predicted to act as effectors, either subverting or triggering host immunity. However, 16 associating particular effectors with virulence or avirulence outcomes presents challenges due to the 17 combinatorial complexity. Here we use defined aphid and host genetics to test for co -segregation of 18 expressed aphid transcripts and proteins with virulent or avirulent phenotypes. 19 Results. We compared virulent and avirulent pea aphid parental genotypes, and their bulk segregant 20 F1 progeny on Medicago truncatula genotypes carrying or lacking the RAP1 resistance quantitative 21 trait locus. Differential gene expression analysis of whole body and head samples, in combination with 22 proteomics of saliva and salivary glands , enabled us to pinpoint proteins associated with 23 virulence/avirulence phenotypes. There was relatively little impact of host genotype, whereas l arge 24 numbers of transcripts and proteins were differentially expressed between parental aphids, likely a 25 reflection of their classification as divergent biotypes within the pea aphid species complex. Many 26 fewer transcripts intersected with the equivalent differential expression patterns in the bulked F1 27 progeny, providing an effective filter for removing genomic background effects. Overall, there were 28 more upregulated genes detected in the F1 avirulent dataset compared with the virulent one. Some 29 genes were differentially expressed both in the transcriptome and in the proteome datasets, with 30 aminopeptidase N prot eins being the most frequent differentially expressed family. In addition, a 31 substantial proportion (27%) of salivary proteins lack annotations, suggesting that many novel 32 functions remain to be discovered. 33 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 2 Conclusions. Especially when combined with tightly controlled genetics of both insect and host, multi-34 omics approaches are powerful tools for revealing and filtering candidate lists down to plausible genes 35 for further functional analysis as putative aphid effectors. 36

Aphid, transcriptomics, proteomics, saliva, effector, virulence, avirulence 37 38

39 Crop losses due to insect pests represent an enduring challenge for agriculture and global food 40 security. Aphids are a major problematic group, due both to the direct damage they cause by phloem 41 sap feeding and to indirect effects through acting as vectors for transmission of many viruses. Impacts 42 of pests are further exacerbated by the breakdown of genetically based crop resistance mechanisms 43 due to selection pressures driving pest evolution, as well as evolved insecticide resistance. 44 In contrast to related fields such as plant -pathogen interactions, the molecular relationships that 45 determine (in)compatibility of plant -aphid interactions are relatively poorly understood. Specific 46 resistance to plant pathogens frequently involves recognition of pathogen effectors, often by 47 resistance proteins (R) characterised by nucleotide -binding and leucine rich repeat (NLR) domains. 48 Several coiled coil domain NLR proteins have been implicated in resistance to aphids and their close 49 relatives. For example, Mi -1, Vat and Bph14 confer resistance to certain biotypes of Macrosiphum 50 euphorbiae (potato aphid) [1], Aphis gossypii (melon-cotton aphid) [2] and Nilaparvata lugens (brown 51 planthopper) [3], respectively. These NLR receptors are predicted to be involved in direct or indirect 52 recognition of molecular signatures that insects, like plant pathogens, release inside their hosts. 53 Indeed, aphids secrete multiple effector proteins into their saliva, that are then predicted to be 54 delivered into plant tissues to modulate host cell processes and to suppress or trigger host defences 55 [4–7]. Although there is one recent report of the BISP effector from brown planthopper, an aphid 56 relative, interacting with the BPH 14 NLR in rice [8], there are currently no examples where cognate 57 aphid effector and NLR pairs have been fully defined. Improved molecular insights into virulence and 58 resistance mechanisms taking place during both compatible and incompatible plant-aphid interactions 59 are therefore a priority, and can provide essential knowledge for future development of durable aphid 60 control strategies. 61 The availability of extensive genome, transcriptome and resequencing resources for the model aphid 62 species Acyrthosiphon pisum (pea aphid) [9, 10] have enabled comprehensive genome -wide 63 explorations. There are also genomic sequences now available at NCBI and Aphid Base 64 (https://bipaa.genouest.org/is/aphidbase/) for more than 25 species of aphids and close relatives , 65 often associated with gene predictions and transcriptomes [11]. In addit ion, several papers have 66 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 3 attempted to define the aphid effectorome, either by direct analysis of salivary proteins, or by 67 transcriptomics of salivary glands, coupled with filters for predicted secreted, non -trans-membrane 68 proteins [12–17]. Beyond the true aphids (superfamily Aphidoidea), there are now genomic resources 69 for sister groups within the Hemiptera such as planthoppers, leafhoppers, psyllids, whitefly and scale 70 insects (https://www.ncbi.nlm.nih.gov/assembly/?term=hemiptera) that likewise are major crop 71 pests, alongside genomes for triatomines and bed bugs, hemipterans that feed on animal rather than 72 plant hosts. Outside the Hemiptera, genomic data have been published for sucking pests such as thrips 73 and spider mites that feed on plant tissues other than phloem [18–20]. Genome, transcriptome and 74 proteome comparisons across clades may enable definition of putative effector subsets that are 75 necessary for different feeding modes, and may provide insights into conserved and divergent modes 76 of action in terms of how the plant immune system is targeted to enable successful parasitism. 77 Despite the wide range of functional genomics studies published to date, one common limitation is 78 the lack of understanding of the differences in effector complements between virulent (host -79 compatible) and avirulent (host -incompatible) genotypes. Genetic differences operate at several 80 taxonomic levels. First, there are major differences across aphid species in their host preferences and 81 host compatibilities. Some species, such as peach potato aphid ( Myzus persicae) are generalists that 82 can feed on at least 40 0 known plant species, making them widespread crop pests [21]. Others are 83 specialists, such as pea aphid (A. pisum) that exclusively feeds on legumes (Fabaceae). Second, there 84 is substantial diversity within species such as A. pisum that has led to its description as a species 85 complex comprising several host races that each have a strong preference for particular legume 86 species, supported by robust molecular marker fingerprints for each host race [22, 23] . T here is 87 evidence of divergence and differential expression of chemosensory gene families such as odorant 88 receptors across different pea aphid biotypes [24, 25] , bu t causative relationships have yet to be 89 established for genes and proteins that govern the range of compatible and incompatible interactions 90 seen. There is also clear evidence that some host races can survive and sometimes thrive as migrants 91 on hosts outside their preferred species range [22]. Finally, at the intra-specific level for both aphids 92 and hosts, there can be a wide range of compatibilities. For example, from testing eight genotypes of 93 A. pisum in combination with 23 different Medicago truncatula (Mt) accessions, we discovered high 94 diversity in both species that did not correspond particularly strongly to host races or to geographic 95 origins of the host lines [26]. Parallel to this, crossing two divergent pea aphid biotypes to generate F1 96 recombinant populations uncovered Mendelian segregation of virulence/avirulence on Mt genotypes 97 carrying the RAP1 aphid resistance QTL [27, 28]. 98 Here, we report global exploration of the molecular basis for aphid virulence and avirulence on 99 defined host genotypes. Specifically, we aimed to link phenotypes to candidate effectors and related 100 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 4 genes by multiple comparisons of the transcriptomes and proteomes of two divergent parental pea 101 aphid clones, along with the transcriptomes of segregating avirulent and virulent pool ed individuals 102 from within F1 cross populations (Fig. 1). We also critically analysed the effectiveness of combined 103 omics approaches as a means to robustly uncover proteins with pivotal biological roles, such as 104 effectors that determine the difference between virulent and avirulent outcomes. 105

and Discussion 106 Generation and analysis of aphid populations for RNA-Seq analyses 107 In our previous work [27], we had demonstrated Mendelian segregation of inheritance of virulent and 108 avirulent phenotypes in F1 pea aphid populations derived from a cross between N116 and PS01 109 (virulent and avirulent parental clones, respectively) when infested on M. truncatula hosts carrying 110 the RAP1 resistance QTL [28]. On this basis, we reasoned that the molecular basis of the difference 111 between virulent and avirulent aphids could be revealed by transcriptomic and proteomic analysis. 112 However, there were likely to be thousands of genetic and gene expression differences between the 113 parental genotypes, that are representatives of phenotypically contrasting biotypes within the highly 114 diverse pea aphid species complex [22, 26] . This makes it difficult to discern unrelated gen omic 115

differences from causative genes responsible for suppressing host immunity or for 116 triggering R -gene dependent defences. To address this challenge, we employed a bulk segregant 117 analysis (BSA-) RNA-Seq approach that would both reduce the genetic background effects and allow 118 us to test for heritability of differentially expressed (DE) genes across parental and F1 generations. 119 Enabling this strategy first required us to re-create the segregating F1 populations previously reported 120 [27]. 121 We induced sexual forms of PS01 and N116 and conducted reciprocal crosses, leading to screening of 122 a total of 78 F1 clones on two host plant genotypes carrying RAP1: Jemalong A17 (hereafter A17), the 123 original source of the identified RAP1 QTL, and a resistant near -isogenic line ( RNIL) derived from a 124 mapping population [29] using A17 as one of the parents. The RAP1 aphid resistance QTL is highly 125 effective against PS01 aphids , typically resulting in high mortality, whereas N116 aphids are 126 unaffected. Progeny were verified as true F1 hybrids by a panel of seven SSR markers [22] and by 127 screening for maternal inheritance of secondary symbionts reported in the pea aphid [30]. Using a 128 virulence index based on a combination of aphid survival and reproduction, F1 clones were first ranked 129 according to performance on A17 . Phenotypes ranged from fully virulent to fully avirulent 130 (Supplementary Material 1A), similar to previous findings [27], although in the present experiment the 131 population as a whole did not display complete segregation into discrete virulent and avirulent 132 categories. As also previously shown, resistance in the RNIL was slightly weaker than in A17, with F1 133 clones ranging from virulent to avirulent, and importantly performance on the two host genotypes 134 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 5 was significantly correlated (Pearson r 0.72, P 1.82 e-13). All F1 clones were virulent on hosts lacking 135 RAP1 (Supplementary Material 1B ). We then selected 22 sibling F1 clones from each end of the 136 distribution to provide two bulk sample sets with the strongest virulent (VIR) and avirulent (AVR) 137 phenotypes for subsequent transcriptomic analysis. Fig. 2 shows the complete separation of the 138 selected clones into virulent and avirulent classifications. As a final check prior to RNA-Seq 139 experiments, we re -confirmed separation of survival rates of these two subsets of clones on both 140 resistant host genotypes (Supplementary Material 1C). 141 Transcriptomic analyses 142 We first ran a n RNA-Seq experiment using the parental clones N116 and PS01 infested onto either 143 A17 or the susceptible DZA315.16 host (hereafter DZA) for 24 h prior to collection of heads for RNA 144 extraction. The multiple aims were to enrich for transcripts from salivary glands that express candidate 145 effectors, to uncover the transcriptome differences between the parental aphid genotypes, and to 146 reveal the impact of host plant genotype. Each aphid x host combination was replicated three times, 147 giving a total of 12 libraries, ranging from 6.8 to 10.6 million reads uniquely mapped to the reference 148 genome (Supplementary Material 2A). 149 Hierarchical c lustering and principal components analysis (PCA) of the transcriptom ic expression 150 profiles both indicated that the replicates of each treatment were closely correlated in all cases, so no 151 datasets needed to be discarded (Fig. 3A,B). These analyses additionally revealed that samples were 152 separated largely by aphid genotype rather than host plant treatment. Overall, the transcriptomes of 153 the two aphid genotypes on A 17 plants were clearly differentiated , with a total of 483 genes 154 significantly upregulated in N116 and 452 in PS01 (log2 fold change >2.0, FDR <0.05; Supplementary 155

3; Fig. 3C). Similarly, on DZA host plants, 395 and 363 genes were upregulated in N116 and 156 PS01, respectively. In contrast, expression of relatively few genes, between three and 27, across all 157 the pairwise comparisons, was significantly affected by the host plant (Supplementary Material 3; Fig. 158 3C). Functions of the DE genes are considered below, in conjunction with the other transcriptomic and 159 proteomic experiments. 160 We next undertook a larger RNA-Seq experiment, sampling whole aphid bodies in order to capture 161 transcripts from all tissues. Using aphids infested onto A17 host plants for 24 h, w e again compared 162 N116 and PS01 parental clones, but this time alongside the bulked segregant pools of VIR and AVR F1 163 clones described above . Five biological replicates for each gave a total of 20 RNA libraries each 164 containing 14 to 22 million reads that uniquely map to the reference genome (Supplementary Material 165 2B). 166 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 6 Similar to the heads experiment, multivariate analysis by hierarchical clustering and PCA both 167 indicated that all replicates within each sample type grouped together, and that each sample type was 168 clearly differentiated. As expected, the genetically divergent parents were again highly separated, 169 whereas the two pooled F1 datasets were much closer to each other , as they contain 50% of each 170 parental genome, with each pool representing the average transcriptome of multiple independent F1 171 clones (Fig. 4A,B). 172 Differentially expressed genes were identified for all pairwise comparisons between samples (Fig. 4C). 173 The number of up and down-regulated genes between the parental pairs and the pair of F1 pools are 174 shown in Fig. 5A, with the gene lists provided in Supplementary Material 3 . Several hundred genes 175 were differentially expressed in both the whole-body and head comparisons of the parents. Some of 176 these DE genes likely reflect genomic differences between the parental clones that are representatives 177 of divergent pea aphid biotypes. However, relatively few DE genes were detected in the F1 samples, 178 with only 24 genes up-regulated in the VIR pool and 64 in the AVR pool. These numbers can also be 179 interpreted as a higher number of genes being down-regulated in the VIR F1 aphids. Fig. 5B,C show 180 the overlaps across head and whole -body datasets for N116 /VIR and PS01 /AVR, respectively . 181 Unexpectedly, the intersections of DE genes revealed subsets where the direction of expression was 182 opposite between the parental pair and the F1 pooled pairs , with three genes upregulated in N116 183 and AVR F1, and 13 genes upregulated in PS01 and VIR F1 (Fig . 5D, Fig. 7G,H). Moreover, very few 184 genes were upregulated in both parental N116 and VIR F1 pool datasets . A plausible explanation is 185 that the genes governing virulence in N116 are not the same as those that result in virulent 186 phenotypes in the F1 population . Each individual in the F1 population carries a random 50% of the 187 genome of each parent , creating a high degree of combinatorial complexit y. Nonetheless, the DE 188 genes in the F1 data derive from the average across the 22 individuals used to create each bulk RNA 189 pool, and are therefore likely to be biologically relevant to virulence or avirulence functions rather 190 than background genomic noise. Such genes merit further exploration in both parental and F 1 191 genotypes. 192 Quantitative proteomic analysis of saliva and salivary glands. 193 To determine whether differences exist between the salivary protein profiles of the two parental 194 aphid clones, a comparative analysis of salivary gland and salivary proteomes was conducted. A total 195 of 2343 and 2276 high confidence proteins were detected from salivary glands of N116 and PS01 , 196 respectively (Supplementary Material 4 ), with 2038 proteins (80%) common to both ( Fig. 6A). Each 197 biotype had similar proportions of non-annotated proteins (PS01: 5.4 % and N116: 6.2%) and proteins 198 predicted to have secretion signals (PS01: 16.6% and N116: 17.3%). These proportions of secreted and 199 non-annotated proteins are typical for pea aphid biotypes [12, 31]. Two major clusters were revealed 200 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 7 by PCA (Fig. 6C), corresponding to the two aphid genotypes. Principal Components 1 and 2 account 201 for 64% of the variation, indicating distinct protein profiles in the salivary glands of each genotype. 202 This distinction was further supported by quantitative analysis that identified 23 5 statistically 203 significant differentially abundant (SSDA) proteins (p<0.05), with 1 36 and 99 proteins having higher 204 abundances in N116 and PS01 salivary glands , respectively ( Fig. 6E; Supplementary Material 4 ). 205 Relative fold changes (RFC) ranged from −48.5 to +140.0 indicating that even when both genotypes 206 engage in compatible interactions with the same plant type ( V. faba in this case) the salivary gland 207 profiles are divergent both qualitatively and quantitatively. 208 Of the 136 SSDA salivary gland proteins with increased abundance in N116, 60 (44%) were predicted 209 to be secreted and 27 (20%) had no annotations. Similar proportions were observed within the 99 210 SSDA proteins with increased abundance in PS01 , with 33 (33%) and 18 (18%) proteins having a 211 secretion signal or no annotations, respectively. The se proportions of secreted and non -annotated 212 proteins within the differentially abundant sets are substantially higher than the corresponding 213 proportions in the background salivary gland proteomes described above. Of the top ten proteins with 214 the highest relative abundance in N116, seven had no annotation : ACPISUM_000319 (ACYPI007553; 215 RFC 140.0) and ACPISUM_029783 ( LOC100573424; RFC 64), ACPISUM_008675 (LOC100162547; RFC 216 32), ACPISUM_016335 (Not annotated; RFC 26), ACPISUM_017388 (LOC103309964; RFC 21.1), 217 ACPISUM_003551 (LOC100534636; RFC, 21.1) and ACPISUM_009099 (LOC112598674, 18.4). The 218 other proteins in the top ten were a kinase ACPISUM_015393 (developmentally -regulated protein 219 kinase 1; RFC 64) and two aminopeptidases (ACPISUM_009259; RFC 36.8 and ACPISUM_005699; RFC 220 22.6). Of the top ten proteins with highest abundances in PS01 in comparison to N116, two were 221 uncharacterised: ACPISUM_007394 (LOC100572241; RFC 48.5) and ACPISUM_007714 222 (LOC100534636; RFC 11.3) ; and two were glutathione S -transferases (ACPISUM_019160 and 223 ACPISUM_001883, both RFCs of 8.6). Other proteins included a different developmentally-regulated 224 protein kinase (ACPISUM_005630; RFC 17.1), a peroxidase (ACPISUM_020816; RFC 9.8), a prostatic 225 spermine-binding protein (ACPISUM_004331; RFC 8), peroxidasin (ACPISUM_019870; RFC 6.5), an 226 ATPase subunit (ACPISUM_009308; RFC 5.7) and glyoxylate reductase (ACPISUM_021751, RFC 4.9). 227 We next examined aphid saliva proteins. Although the samples are collected from artificial diets, these 228 salivary secretomes are likely to be highly similar to the proteins delivered into plant tissues during 229 interactions with the host, and therefore are predicted to include the entire set of effectors. We 230 focussed on categorisation of the total salivary protein lists, and of the DE proteins. Although the 231 analysis of saliva revealed far fewer proteins than from the salivary gland samples, there is again a 232 clear distinction between the two genotypes. A total of 69 and 97 high confidence proteins were found 233 in N116 and PS01 saliva, respectively (Fig. 6B; Supplementary Material 4) with 22 (32% for N116) and 234 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 8 50 (52% for PS01) proteins, being deemed unique to each. A large proportion (30% for PS01 and 25% 235 for N116) of the salivary proteomes had no annotations, indicating their potential phylogenetic 236 restriction to aphids. In addition, 39% and 3 2% of the proteins had predicted canonical secretion 237 signals for PS01 and N116 saliva, respectively. Notably, although s aliva proteins detected in diet 238 samples have, by definition , been secreted , the majority appear not to have canonical secretion 239 signals. Explanations range from incomplete/incorrect gene models to non-canonical or alternative 240 secretion mechanisms. Our results highlight the importance of combining several approaches when 241 attempting to identify potential effectors and molecular determinants of virulence/avirulence. 242 Omitting proteins without secretion signals from bioinformatic pipelines may result in many effector 243 candidates being overlooked. 244 As with the salivary glands, PCA of the salivary proteins completely resolved two groups, with PC1 and 245 PC2 accounting for 94% of the total variation (Fig. 6D). Label free quantitative analysis using MaxQuant 246 identified 47 SSDA proteins with 12 and 35 proteins having higher abundance in N116 and PS01 saliva, 247 respectively ( Fig. 6F; Supplementary Material 4 ). Notably, N116 saliva comprises fewer detected 248 proteins and fewer SSDA proteins than PS01, possibly pointing to a strategy that enables evasion of 249 host defences. If , for example, one or more of the proteins uniquely detected in PS01 saliva act as 250 avirulence factors due to cognate receptors in the host plant, their absence or low abundance in N116 251 may result in a compatible interaction. However, it remains to be experimentally determined whether 252 these genotypic differences in type or number of saliva proteins are causatively associated with 253 virulence or avirulence. 254 Most of the salivary proteins identified here have previously been associated with pea aphid saliva 255 including multiple members of M1 and M2 metalloprotease families , along with peroxidases, 256 glutathione-S-transferases, glucose dehydrogenase and regucalcin [12, 32] . Apart from the 257 Aminopeptidase N (APN) category discussed in detail below, the most frequent annotation was for 258 unknown proteins: 20-26% of the total saliva list for each clone, and 21% of the DE saliva proteins. 259 Four out of the ten DE unknown proteins also featured within the top 20 proteins by MS intensity or 260 protein coverage. High proportions of unknown proteins have been noted in earlier studies of aphid 261 saliva and the salivary gland predicted secretome [31]. In addition, a homologue of a salivary effector 262 previously characterised for Myzus persicae (Mp1) [33] had a higher abundance in PS01 saliva 263 (ACPISUM_000421; RFC 14). The relative fold changes of salivary proteins ranged from -2352 for 264 regucalcin to 724 for members of the APN (M1 metalloprotease) family, which represented the most 265 differentially abundant proteins in PS01 and N116 saliva, respectively. Although these RFC values can 266 be considered arbitrary due to imputation of low abundant values in samples where the proteins are 267 in fact absent, the re is very clear divergence of salivary proteomes both in the proteins uniquely 268 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 9 detected in one or other genotype , and in the large differences in apparent abundance of several 269 proteins present in both genotypes. The full lists of proteins exclusively found in the saliva or salivary 270 gland proteomes of both genotypes are provided in Suppleme ntary Material 4 , with 25 and five 271 proteins exclusive to the salivary glands and saliva of N116, respectively. For PS01, the corresponding 272 numbers were 10 and 13 proteins exclusive to the salivary glands and saliva, respectively. These 273 proteins were present in all replicates of one genotype while being absent in all replicates of the other, 274 strongly supporting their status as candidate effectors, that may individually or collectively determine 275 the VIR and AVR phenotypes observed for each genotype on different host plants. 276 Comparison of the quantitative differences in protein abundance across both the saliva and salivary 277 gland datasets revealed clear similarities in the two proteomes analysed for each genotype . Five 278 proteins that were of higher abundance in N116 saliva were also more abundant in N116 salivary 279 glands in comparison to their PS01 counterparts. A similar trend was observed for nine PS01 salivary 280 and salivary gland proteins ( Supplementary Material 4 ), with the RFCs for these proteins positively 281 correlated across both biological sample types. The fact that the abundances of these salivary gland 282 proteins are mirrored at the level of externally delivered oral secretions highlights the robustness of 283 both analyses, and points to likely roles as virulence or avirulence determinants in two genotypes with 284 distinct host preferences. Such proteins represent excellent candidates for future characterisation to 285 determine their effector status , especially those that are also supported by DE transcript profiles 286 (Table 1). 287 Overlap between transcriptomics and proteomics datasets 288 Across the transcriptomics and prote omics experiments, we analysed all the intersections then 289 extracted the proteins and DE gene subsets that showed the greatest overlaps (Table 1; 290 Supplementary Material 3 and 4), partitioning into genes/proteins associated with virulence, in N116 291 or the VIR F1 pool, or with avirulence, in PS01 or the AVR F1 pool. The number of DE genes or proteins 292 in the hea d transcriptome, whole body transcriptome and salivary gland proteome datasets were 293 broadly similar between VIR and AVR samples. However, the PS01 saliva protein and the AVR F1 pool 294 transcript lists were longer than those for N116 saliva and VIR F1 pool transcripts, reflected by larger 295 intersections in the former. Over half (33/64) of genes upregulated in the AVR F1 pool were also in at 296 least one other list, whereas only three out of 24 intersected from the VIR F1 pool data. Whole body 297 RNA-Seq data for a selection of these intersected genes are plotted in Fig . 7. Several of the AVR -298 upregulated genes shown are annotated as enzymes with hydrolase, glycosidase or peroxidase 299 functions. Other annotations include a transcription factor and proteins of unknown function. Genes 300 on the VIR side included ACPISUM_013796 (myrosinase 1 -like) and ACPISUM_019971 (glutathione 301 hydrolase 1 proenzyme -like), although these were not found in saliva . Across the multiple 302 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 10 experiments, the two most frequently found genes in the AVR data were ACPISUM_021997 303 (regucalcin-like) previously reported as a Ca-binding protein [32], present in all lists except heads RNA, 304 and ACPISUM_029930 (uncharacterized protein LOC100575698), present in all five lists. These AVR-305 related salivary proteins represent strong candidates for functional effectors, based on the multiple 306 strands of evidence for their differential expression and importantly for co-segregation of their 307 expression with the avirulence phenotype in the F1 population . We have therefore uncovered 308 heritable differences in salivary proteins that associate with avirulence, in this case an incompatible 309 phenotype on Mt hosts carrying the RAP1 QTL [27, 28]. Intriguingly, however, we found no equivalent 310 strong candidates for salivary proteins that might represent the dominant virulence factor predicted 311 by previous genetic studies [27]. Alternative explanations for the Mendelian segregation found in that 312 study could be that the proposed “virulence” gene is not an effector per se, but instead could be an 313 upstream positive regulator, or a negative regulator of one or more effectors that act as avirulence 314 factors detected by a RAP1 dependent pathway. 315 Gene Ontology analysis 316 We undertook Gene Ontology (GO) analysis to reveal functional categories and genes that were 317 enriched in the differentially expressed gene and protein data sets. Using a FDR of <0.05, many gene 318 sets contained few or no significant ly enriched terms (Table 2; Supplementary Material 5). For the 319 whole-body transcriptome data, aminopeptidase N (APN) proteins were strongly enriched, with 320 different genes within this family upregulated in each of the parental aphids (discussed further below). 321 These trends were reinforced by comparison of parental transcriptomes in the heads RNA-Seq 322 analyses where APN proteins were similarly enriched in both parents. The DE gene sets between the 323 pooled VIR and AVR F1 samples indicated no enriched terms in the VIR data, and only a single term 324 among the AVR upregulated genes: glucosidase II complex, that localises to the ER. These two gene 325 sets are both relatively small (64 and 24 genes), reducing the likelihood of finding significant trends. 326 Because very few significantly enriched terms were revealed by the initial GO analyses, we applied a 327 lower stringency to inform wider trends in each of the DE gene sets. Here, we examined all terms for 328 which at least two genes and a significant P value (<0.05) were returned. For the DE gene sets from 329 RNA-Seq of heads, the majority of enriched terms were associated with the virulent N116 parent on 330 both host genotypes. Although there was obvious redundancy of many terms, a substantial proportion 331 (30-40%) for N116 relate to energy metabolism including mitochondria, TCA cycle, oxidative 332 phosphorylation and lipid metabolism. In contrast, the PS01 enriched terms included several for 333 protein processing including peptidases, proteolysis and protein glycosylation ; and several for ATP-334 related transport (Supplementary Material 5 ). When each parental aphid genotype was compared 335 separately for its differential responses to the two host genotypes (A17 and DZA), no significant terms 336 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 11 were found for PS01, and only one weakly significant term for N116: polytene chromosome puffing. 337 The equivalent GO analysis of whole body RNA-Seq data returned significantly enriched terms for both 338 aphid genotypes, including several for protein modification (Supplementary Material 5). 339 For the DE datasets from salivary gland proteomes, the lower stringency analysis revealed enrichment 340 of distinct functional categories for each parental genotype. For N116, protein modification terms 341 were prevalent including peptidase activity, serine -type endopeptidase inhibitor activity, negative 342 regulation of protein metabolic process, aminopeptidase activity, protein kinase binding and 343 regulation of protein phosphorylation. In contrast, for PS01, ATPase terms were predominant 344 including several related to membrane transport , as also found in the PS01 heads RNA -Seq data 345 (Supplementary Material 5). 346 Exopeptidases are abundant in saliva, and the majority are DE between aphid genotypes 347 The saliva protein total and DE lists were much shorter, precluding formal GO analysis, but manual 348 inspection indicated high proportions of exopeptidases: a total of 29 different proteins (Table 3) , 349 representing 22-34% of the protein list for each genotype. These were mainly APN proteins but also 350 four members of the Angiotensin Converting Enzyme (ACE) family that are M2 metalloproteases with 351 carboxypeptidase activity. The abundance of APNs in the saliva protein list broadly corroborates the 352 major enriched GO categories detected in the transcriptome analyses. 353 Most of the exopeptidases detected from aphid saliva ( 23/29; 7 9%) were differentially abundant 354 between the parental aphid genotypes. Twenty-two of the 29 saliva exopeptidases were also found in 355 the salivary gland proteomes, with many showing the same direction of differential expression (9 APN, 356 2 ACE). Moreover, 1 5 (60%) of th e APN proteins were DE in heads and /or whole body RNA-Seq 357 samples (Table 3). Previous reports on pea aphid saliva and salivary gland components have also 358 reported multiple APN and ACE proteins [12, 13, 32, 34]. Similar to our findings, one of these studies 359 reported 11 APN genes that were differentially expressed in a biotype -specific manner, with five of 360 these detected as proteins in saliva [13]. Taking all the evidence together, it is clear that the APN family 361 is highly diversified in pea aphids and represents a major component of the salivary proteome by 362 several measures: the high total number of proteins detected , many of these proteins are high 363 abundance (13 of 20 top scoring in both N116 and PS01 saliva), and most are differentially expressed 364 between aphid genotypes. 365 Aphid and mammalian ACE proteins have similar sequences and may have broadly similar functions 366 as dipeptidases or by cleaving a single amino acid from the C terminus. However, mammalian ACE 367 proteins are membrane anchored whereas aphid ACEs carry secretion signals, consistent with their 368 detection in saliva. The exact catalytic functions and biological roles of aphid ACE and APN proteins 369 remain to be determined. Cleavage of proteins and peptides could relate to targeting host proteins 370 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 12 such as those involved in defensive sieve -tube blocking as shown at least for the atypical 371 extrafascicular phloem exudate of cucurbits [35]. Alternatively, although there is currently no direct 372 evidence, exopeptidases may act on other salivary protein components, for example to process 373 effectors into active forms. Another non-mutually exclusive possibility is a role in aphid nutrition, with 374 many insects using extra-organismal (extra-oral) digestion typical of arthropods including Hemiptera, 375 enabling nutrition capture from large host s/prey [32, 36] . Exopeptidases typically release N or C 376 terminal single amino acids and dipeptides, potentially enabling supply of essential amino acids, some 377 of which cannot be biosynthesised directly from the enzyme repertoires of hemimetabolous aphids. 378 Multi-omic approaches to detecting candidate effectors 379 We compared the efficiencies of the four different experiments in terms of detecting aphid candidate 380 effectors and related genes: RNA-Seq of heads and whole bodies, and proteomics of saliva and salivary 381 glands. For all datasets, we focussed mainly on differential expression between the highly divergent 382 parental clones N116 and PS01. Because saliva represents the “ground truth” of proteins predicted to 383 be delivered into plant host tissues, we additionally considered saliva proteins that were detected but 384 not DE. Although the proteomics methods are highly sensitive, there are likely to be some further low 385 abundance salivary proteins that were not detected here. In addition, there may be some salivary 386 proteins that are only expressed in response to aphids interacting with their host plants , and hence 387 would not be found in artificial diet samples. Likewise , some proteins may not be stable under the 388 artificial diet conditions. As a case study, we selected the significantly enriched exopeptidases , that 389 comprised the large APN family and the smaller group of ACE proteins. We compared success of 390 detecting genes from the saliva data in the other three experiments, and noted whether the same DE 391 patterns were found (Table 3). The overall trends were broadly correlated, with 18/24 (75%) DE saliva 392 proteins also found to be DE in at least one of the other approaches. Only two genes showed a 393 mismatch in DE direction: ACPISUM_009259 between salivary gland and whole body; and 394 ACPISUM_020790 between saliva and salivary gland. Individually, RNA-Seq of heads was the most 395 effective experiment (14/24) at corroborating the DE saliva protein evidence, followed by RNA-Seq of 396 whole bodies (10) and proteomics of salivary glands (8). 397 There are several reports where effectors are predicted f rom aphid salivary gland transcriptomes or 398 proteomes, or other transcriptome datasets, typically filtering for presence of a signal peptide or other 399 secretion motif, and absence of transmembrane domains [12–17]. For our exopeptidase data (Table 400 3), we detected an additional seven APNs in salivary gland proteomes or the transcriptome data, that 401 were not found in saliva, of which five were DE in at least one dataset. The ir absence from saliva 402 indicates these proteins may be considered false positives for candidate effectors, although some low 403 expressed proteins may go undetected . We considered which of the approaches was the most 404 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 13 effective at detecting candidate effectors, and whether multiple omics approaches are advantageous, 405 noting that all require substantial resource investment. Although saliva collection is an exacting and 406 time-intensive procedure, saliva proteomics provided the greatest coverage of candidate effectors 407 here, and quantitative analysis of mass spectrometry data enables robust assignment of differential 408 expression. Of the other approaches, RNA-Seq of heads may be the most effective means to 409 complement the saliva analyses by reinforcing evidence of differential expression, but in the work 410 here did not greatly extend the effector lists per se. 411

412 In this study, we demonstrated that transcriptomics and proteomics are both highly effective tools for 413 discovering differentially expressed aphid genes and proteins . The protein subsets present in saliva 414 are likely candidates for effectors with virulence and/or avirulence functions in host plants , and 415 represent priorities for further study especially to determine if differential protein abundance is 416 inherited into the segregating F1 aphid populations. Precise biochemical functions and host targets of 417 most of these effectors are also currently unknown even in cases, such as the exopeptidases, where 418 there are confident gene annotations. Exopeptidases are dominant in saliva by number of different 419 proteins, by frequency of differential abundance, and by quantity. Likewise, there are many proteins 420 of unknown function, with a substantial proportion found at high levels in saliva. Some of these 421 unknown proteins may prove to be pivotal in explaining aphids’ unique and highly successful lifestyle 422 as phloem feeders. 423

424 Aphids and crossing 425 Pea aphid (Acyrthosiphon pisum) clones were maintained on tic bean (Vicia faba minor) as described 426 in [26]. Parental genotypes were PS01 and N116. PS01 is a biotype adapted to Pisum sativum whereas 427 N116 is adapted to Medicago sativa [26]. Reciprocal crosses were made between PS01 and N116 to 428 generate F1 hybrid populations, following the protocol of [27]. In brief, parthenogenetic females were 429 induced to generate sexual forms by transfer to short days and lower temperatures to simulate 430 autumn. Eggs resulting from controlled matings were collected onto moist filter paper in petri dishes, 431 and subjected to 90 to 105 days at 4°C to induce exit from diapause. Individual hatchlings were 432 subsequently used to generate multiple parallel clonal F1 lineages. Parents and progeny were 433 genotyped with a set of seven microsatellite markers [22] to verify correctness of crosses. All new F1 434 progeny were maintained for at least three generations before testing performance on different host 435 plants. 436 Plants and assessment of virulence 437 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 14 Based on previous findings [27], PS01 aphids are avirulent on Medicago truncatula J A17 that carries 438 the resistance QTL, RAP1 [28]. Near isogenic lines (NILs) derived from a cross (LR4 [29]) between A17 439 and M. truncatula DZA315.16 were also used. PS01 is likewise incompatible with the resistant NIL 440 (RNIL), but is compatible with the susceptible NIL (SNIL) and with DZA315.16. N116 aphids are 441 compatible with all these genotypes. F1 progeny were tested for virulence on both A17 and RNIL, 442 based on [26]. Briefly, five nymphs of each clone were infested onto ten A17 or RNIL plants, then 443 scored for survival and production of new nymphs 10 d later. At least 40 F1 clones each of PS01 x N116 444 and N116 x PS01 were screened. An overall virulence index was adapted from a calculation proposed 445 in [37]: 446 Virulence index = log2 (mean number surviving out of 5 x number of nymphs produced + 1) 447 Virulent (VIR) clones were defined as index >4 and >5 on A17 and RNIL, respectively, and avirulent 448 (AVR) clones were correspondingly defined as index <2 and <4. The different category thresholds on 449 A17 and RNIL reflect the latter’s slightly lower resistance. Clones falling into the same phenotype 450 category (VIR or AVR) on both A17 and RNIL were then subject to a further confirmation screen where 451 survival on A17 and RNIL was counted 5 d after infestation. In the confirmation experiment, four plants 452 were used for each aphid x host combination, with five aphids infested onto each plant. Cutoffs were 453 >80% survival for virulence on both hosts, and <20% and <70% for avirulence on A17 and RNIL, 454 respectively. A few F1 clones showed relatively high survival at 5 days but had very weak growth, and 455 therefore were categorised as AVR. Only F1 clones displaying the same phenotype category on all 456 screening experiments were used subsequently in molecular experiments. 457 Sampling for RNA-Seq 458 Heads experiment: Young adult aphids of clones N116 and PS01, cultured on Vicia faba minor, were 459 infested onto either A17 or DZA315.16 M. truncatula plants for 24 h, then heads (40 per sample) were 460 dissected and frozen immediately on dry ice then stored at -80°C. Three replicates were done for each 461 aphid x plant combination. 462 Whole body experiment: Samples were parental aphid clones (N116 and PS01) and pools of VIR and 463 AVR F1 progeny. Aphids of each individual genotype, age 2 to 3 d, were placed on independent A17 464 plants for 24 h then frozen in liquid nitrogen and stored at -80°C until processing. A total of 22 VIR and 465 22 AVR F1 aphid clones were collected individually, before pooling five aphids of each genotype to 466 comprise one sample. Five biological replicates were analysed for both parental and pooled F 1 467 genotypes. 468 RNA extraction, library preparation and sequencing 469 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 15 Heads were dissected and processed as described in [16]. Total RNA was extracted using a plant RNA 470 extraction kit (Sigma -Aldrich). Illumina TruSeq stranded mRNA -Seq libraries were sequenced at the 471 Genome Sequencing Unit at the University of Dundee on an Illumina HiSeq 2000. 472 RNA for the BSA-RNA-Seq analysis was isolated from three two to three day old nymphs of parental 473 lines (N116, PS01), 22 VIR F1 lines and 2 2 AVR F1 lines, using the Norgen Plant and Fungal RNA kit 474 (Sigma E4913). The RNA isolation followed the instructions of the company supplementing Lysis buffer 475 C with ß-mercaptoethanol. An on -column DNase digest was performed (RNase -Free DNase Set, 476 Qiagen) and the concentration of each sample determined via a Qubit fluorometer with the QubitTM 477 RNA Broadrange (BR) assay kit (Thermo Fisher Scientific). Samples corresponding to five replicates of 478 each of the parental lines and the VIR and AVR F1 pools were used to generate a total of 20 Illumina 479 TruSeq stranded mRNA-Seq libraries which were sequenced in 150 bp paired-end mode on an Illumina 480 HiSeq4000 at Edinburgh Genomics. 481 RNA-Seq data processing and visualisation. 482 Illumina RNA sequence reads were subjected to quality control using FastQC. The reads were the 483 trimmed using Trimmomatic (version 0.32) Q15, min length 55. The trimmed fastq files were the n 484 quasi mapped to the nucleotide gene sequences for the pea aphid using salmon version 1.1. For the 485 pilot study, STAR (2.4.1b) [38] was used to map the reads to the pea aphid genome and HTseq counts 486 was used to quantify the gene expression using AphidBase_OGS2.1b gene annotations. 487 Clone-specific de novo RNA-Seq assemblies (from both the heads and whole -body studies) were 488 individually and collectively generated using Trinity version 2.9.1. All the data were pooled into one 489 for the “collective” assembly, which was used for transcript differential expression analysis. The 490 individual assemblies were used for gene prediction at a later stage. All RNA -Seq assemblies were 491 quality filtered using Transrate to reduce the probability of mis -assembled transcripts. Predicted 492 coding sequences were generated using TransDecoder (with PFAM and BLAST guides). Diamond was 493 used to search against GenbankNR database. Differential expression analysis was performed using 494 EdgeR. Heatmaps and expression profile clustering w ere generated using the ptr script from within 495 the Trinity package. 496 During early analysis, following visual assessment of RNA -seq read mapping and initial differential 497 expression results, we found that the original pea aphid gene predictions (AphidBase_OGS2.1b) and 498 the gene predictions fro m [39] did not fully match those generated by the de novo transcriptome 499 assemblies. Therefore, gene annotation was re -predicted on the published pea aphid genome 500 (OGS2.1b) to improve the accuracy of the gene models. Funannotate, in Other Eukaryotic mode, was 501 used to predict the genes using the de novo RNA-Seq assembly generated above, with RNA-Seq data 502 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 16 mapped using STAR (see above). A total of 29,930 genes were assigned codes in the format 503 ACPISUM_0xxxxx, with the annotations provided at doi.org/10.5281/zenodo.11103500 [40]. 504 To assign the various gene call s from the original genome assembly, bedtools intercept was used to 505 identify genes with overlapping coordinates. If the genes overlapped, then they were considered the 506 same gene. A simple BLAST appro ach could not be used here due to the duplicated nature of aphid 507 assemblies. A combination of reciprocal best BLAST hit, Orthofinder and MCL clustering were used to 508 assign genes between the clones as orthologues. 509 Saliva Collection 510 For proteomics samples, N116 and PS01 were maintained separately on Vicia faba c.v. The Sutton, 511 grown in standard potting compost and kept at 20 oC and a photoperiod of 16 -h light/8 -h dark. 512 Approximately 3,000 mixed aged aphids were positioned on 30 perspex rings (radius 4.5 cm, height 5 513 cm), each containing 4.5 ml of a chemically-defined diet, formulation A from [41], held between two 514 stretched sheets of ParafilmTM. The aphids were reared on the diets at 20°C with 18h light and 6h dark 515 for 24 h after which the diets were pooled and collected and stored at -80°C until required. Four 516 independent replicates were produced by pooling the collected diet from two daily collections 517 (approximately 150 ml). Pooled diets were concentrated using a Vivacell 250 Pressure Concentrator 518 (Sartorius Mechatronics, UK) using a 5000 Da molecular weight cut -off (MWCO) polyethersulfone 519 (PES) membrane. When the final volume had reached 5 ml it was removed and 1 ml of filtered 520 sterilised PBS (phosphate-buffered saline) supplemented with Roche cOmplete TM protease inhibitor 521 cocktail (PIC) was added. The resulting mixture was further concentrated to approximately 250 μl 522 using a Vivaspin 6 centrifuge concentrator (Sartorius Mechatronics, UK) with a 5000 Da MWCO PES 523 membrane, purified using a 2D Clean -up Kit (GE HealthCare) following the manufacturer’s 524 instructions. The resulting protein pellet was suspended in 25 μl 6 M urea, 2 M thiourea, 0.1 M Tris-525 HCl, pH 8.0 and re-quantified using the Qubit Fluorometer. Four independent biological replicates per 526 genotype were subjected to mass spectrometry. 527 Salivary glands 528 The salivary glands from 14-16 day old adult aphids of N116 and PS01 were dissected in ice-cold PBS 529 and transferred to 60 µl PBS supplemented with PIC. Forty pairs of salivary glands were pooled per 530 replicate and homogenized with a motorised, disposable pestle. Sixty microliters of 12 M urea, 4 M 531 thiourea, and PIC was added and the samples were homogenised further and centrifuged at 9,000 × g 532 for 5 min to pellet cellular debris. The supernatant was removed and quantified, and 100 µg of protein 533 was purified using a 2D Clean -up Kit (GE HealthCare) following the manufacturer’s instructions with 534 the exception that 400 μl of precipitant and co-precipitant were used in the first step . The resulting 535 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 17 protein pellet was re-suspended in 30 μl 6 M urea, 2 M thiourea, 0.1 M Tris -HCl, pH 8.0 and re -536 quantified using the Qubit Fluorometer. Four biological replicate s per genotype were subjected to 537 mass spectrometry. 538 Protein sample digestion for mass spectrometry 539 The digestion protocol was the same for both saliva and salivary gland samples and involved the 540 addition of 50 μl ammonium bicarbonate, reduction with 0.5 M dithiothreitol at 56°C for 20 min and 541 alkylation with 0.55 M iodoacetamide at room temperature for 15 min, in the dark. One μl of a 1% 542 w/v solution of ProteaseMax Surfactant Trypsin Enhancer (Promega) and 1 μg of Sequence Grade 543 Trypsin (Promega) were added , then samples were incubated at 37°C for 18 h. Digestion was 544 terminated by adding 1 μl of 100% trichloroacetic acid (Sigma Aldrich) and incubati ng at room 545 temperature for 5 min. Samples were centrifuged for 10 min at 13,000 x g and the supernatant was 546 removed to new microcentrifuge tubes. 547 Mass spectrometry and proteomic data analysis 548 One μg of digested peptide was loaded onto a Dionex Ultimate 3000 (RSLCnano) chromatography 549 system connected to a QExactive (ThermoFisher Scientific) high -resolution accurate mass 550 spectrometer. Peptides were separated by an increasing acetonitrile gradient on a Biobasic C18 551 PicofritTM column (100 mm length, 75 µm ID), using 120 and 50 min reverse phase gradient s for 552 salivary glands and saliva, respectively, at a flow rate of 250 nl min-1. All data were acquired with the 553 mass spectrometer operating in automatic data dependent switching mode. A high -resolution MS 554 scan (300 -2000 Da) was performed using the Orbitrap to select the 15 most intense ions prior to 555 MS/MS. 556 Protein identification and normalisation was conducted using the Andromeda search engine in 557 MaxQuant (version 1.6.17.0; http://maxquant.org/) to correlate the data against the predicted 558 protein set generated in this study (ACPISUM_Proteins; 30891 entries) using default search 559 parameters for Orbitrap data. False Discovery Rates were set to 1% for both peptides and proteins 560 and the FDR was estimated following searches against a target -decoy database. Two searches were 561 conducted for both N116 and PS01 saliva and salivary glands. The first involved a combined search of 562 the raw files for each genotype separately to generate comprehensive proteomes for the saliva or 563 salivary gland (hereafter All Identified Proteins). The second involved a quantitative search of the raw 564 files for all biological replicates (n=4) for the saliva or salivary glands. Quantitative and statistical 565 analyses were conducted in Perseus (Version 1.6.1.1 http://maxquant.org/) using the n ormalized 566 label-free quantitation ( LFQ) intensity values from each sample . The data were filtered to remove 567 contaminants, and peptides identified by site. LFQ intensity values were log2 transformed and samples 568 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 18 were allocated to their corresponding gr oups. A data imputation step was conducted to replace 569 missing values with values that simulate signals of low abundant proteins chosen randomly from a 570 distribution specified by a downshift of 2.1 times the mean standard deviation (SD) of all measured 571 values and a width of 0.1 times this SD. Normalized intensity values were used for principal 572 components analysis. A two-sample t-test was performed using a cut-off value of p ≤ 0.05 to identify 573 statistically significant differentially abundant (SSDA) proteins. Volcano plots were produced by 574 plotting –Log p-values on the y-axis and Log2 fold-change values on the x-axis to visualize differences 575 in protein abundance between the two genotypes. 576 Gene annotations and Gene Ontology analysis 577 Secretion signal properties were predicted using SignalP4.1 [42]. Non-annotated genes were defined 578 as those with the following descriptors : hypothetical protein, uncharacterized protein , NA or 579 ACYPIxxxxx without any other assigned function . GO enrichment analyses w ere performed using 580 GOseq [43]. 581 Data availability 582 Genome annotations: zenodo.org/records/11103500 [40] 583 RNA-Seq: Pea aphid clones N116 and PS01 reared on Medicago truncatula A17 and DZA 315.16, 584 dissected heads: BioProject PRJNA757589, ncbi.nlm.nih.gov/bioproject/PRJNA757589/ 585 RNA-Seq: Pea aphid clones N116 , PS01 and bulk F1 hybrid progeny reared on Medicago truncatula 586 A17, whole body samples: BioProject PRJNA757896, ncbi.nlm.nih.gov/bioproject/PRJNA757896 587 Scripts: github.com/peterthorpe5/Pea_aphid_on_medicago_DZA_A17 588 Proteomics: mass spectrometry data have been deposited to the ProteomeXchange Consortium via 589 the PRIDE partner repository [44], dataset identifiers PXD053355 and PXD053620. 590 591 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 19

592 1. Rossi M, Goggin FL, Milligan SB, Kaloshian I, Ullman DE, Williamson VM. The nematode resistance 593 gene Mi of tomato confers resistance against the potato aphid. Proc Natl Acad Sci. 1998;95:9750–4. 594 2. Dogimont C, Chovelon V, Pauquet J, Boualem A, Bendahmane A. The V at locus encodes for a CC ‐ 595 NBS ‐ LRR protein that confers resistance to A phis gossypii infestation and A . gossypii ‐mediated 596 virus resistance. Plant J. 2014;80:993–1004. 597 3. Du B, Zhang W, Liu B, Hu J, Wei Z, Shi Z, et al. Identification and characterization of Bph14 , a gene 598 conferring resistance to brown planthopper in rice. Proc Natl Acad Sci. 2009;106:22163–8. 599 4. Rodriguez PA, Bos JIB. Toward Understanding the Role of Aphid Effectors in Plant Infestation. Mol 600 Plant-Microbe Interactions®. 2013;26:25–30. 601 5. Elzinga DA, De Vos M, Jander G. Suppression of Plant Defenses by a Myzus persicae (Green Peach 602 Aphid) Salivary Effector Protein. Mol Plant-Microbe Interactions®. 2014;27:747–56. 603 6. Wang W, Dai H, Zhang Y, Chandrasekar R, Luo L, Hiromasa Y, et al. Armet is an effector protein 604 mediating aphid‐plant interactions. FASEB J. 2015;29:2032–45. 605 7. Mugford ST, Barclay E, Drurey C, Findlay KC, Hogenhout SA. An Immuno-Suppressive Aphid Saliva 606 Protein Is Delivered into the Cytosol of Plant Mesophyll Cells During Feeding. Mol Plant-Microbe 607 Interactions®. 2016;29:854–61. 608 8. Guo J, Wang H, Guan W, Guo Q, Wang J, Yang J, et al. A tripartite rheostat controls self-regulated 609 host plant resistance to insects. Nature. 2023;618:799–807. 610 9. The International Aphid Genomics Consortium. Genome Sequence of the Pea Aphid 611 Acyrthosiphon pisum. PLoS Biol. 2010;8:e1000313. 612 10. Li Y, Park H, Smith TE, Moran NA. Gene Family Evolution in the Pea Aphid Based on 613 Chromosome-Level Genome Assembly. Mol Biol Evol. 2019;36:2143–56. 614 11. Shigenobu S, Yorimoto S. Aphid hologenomics: current status and future challenges. Curr Opin 615 Insect Sci. 2022;50:100882. 616 12. Boulain H, Legeai F, Guy E, Morlière S, Douglas NE, Oh J, et al. Fast Evolution and Lineage-Specific 617 Gene Family Expansions of Aphid Salivary Effectors Driven by Interactions with Host-Plants. Genome 618 Biol Evol. 2018;10:1554–72. 619 13. Boulain H, Legeai F, Jaquiéry J, Guy E, Morlière S, Simon J-C, et al. Differential Expression of 620 Candidate Salivary Effector Genes in Pea Aphid Biotypes With Distinct Host Plant Specificity. Front 621 Plant Sci. 2019;10:1301. 622 14. Wang D, Yang Q, Hu X, Liu B, Wang Y. A Method for Identification of Biotype-Specific Salivary 623 Effector Candidates of Aphid. Insects. 2023;14:760. 624 15. Bos JIB, Prince D, Pitino M, Maffei ME, Win J, Hogenhout SA. A Functional Genomics Approach 625 Identifies Candidate Effectors from the Aphid Species Myzus persicae (Green Peach Aphid). PLOS 626 Genet. 2010;6:e1001216. 627 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 20 16. Thorpe P, Cock PJA, Bos J. Comparative transcriptomics and proteomics of three different aphid 628 species identifies core and diverse effector sets. BMC Genomics. 2016;17:172. 629 17. Nicolis VF, Burger NFV, Botha A-M. Whole-body transcriptome mining for candidate effectors 630 from Diuraphis noxia. BMC Genomics. 2022;23:493. 631 18. Rotenberg D, Baumann AA, Ben-Mahmoud S, Christiaens O, Dermauw W, Ioannidis P, et al. 632 Genome-enabled insights into the biology of thrips as crop pests. BMC Biol. 2020;18:142. 633 19. Grbić M, Van Leeuwen T, Clark RM, Rombauts S, Rouzé P, Grbić V, et al. The genome of 634 Tetranychus urticae reveals herbivorous pest adaptations. Nature. 2011;479:487–92. 635 20. Huo S-M, Yan Z-C, Zhang F, Chen L, Sun J-T, Hoffmann AA, et al. Comparative genome and 636 transcriptome analyses reveal innate differences in response to host plants by two color forms of the 637 two-spotted spider mite Tetranychus urticae. BMC Genomics. 2021;22:569. 638 21. Mathers TC, Chen Y, Kaithakottil G, Legeai F, Mugford ST, Baa-Puyoulet P, et al. Rapid 639 transcriptional plasticity of duplicated gene clusters enables a clonally reproducing aphid to colonise 640 diverse plant species. Genome Biol. 2017;18:27. 641 22. Peccoud J, Ollivier A, Plantegenest M, Simon J-C. A continuum of genetic divergence from 642 sympatric host races to species in the pea aphid complex. Proc Natl Acad Sci. 2009;106:7495–500. 643 23. Caillaud MC, Mondor‐Genson G, Levine‐Wilkinson S, Mieuzet L, Frantz A, Simon JC, et al. 644 Microsatellite DNA markers for the pea aphid Acyrthosiphon pisum. Mol Ecol Notes. 2004;4:446–8. 645 24. Eyres I, Jaquiéry J, Sugio A, Duvaux L, Gharbi K, Zhou J, et al. Differential gene expression 646 according to race and host plant in the pea aphid. Mol Ecol. 2016;25:4197–215. 647 25. Eyres I, Duvaux L, Gharbi K, Tucker R, Hopkins D, Simon J, et al. Targeted re‐sequencing confirms 648 the importance of chemosensory genes in aphid host race differentiation. Mol Ecol. 2017;26:43–58. 649 26. Kanvil S, Powell G, Turnbull C. Pea aphid biotype performance on diverse Medicago host 650 genotypes indicates highly specific virulence and resistance functions. Bull Entomol Res. 651 2014;104:689–701. 652 27. Kanvil S, Collins CM, Powell G, Turnbull CGN. Cryptic Virulence and Avirulence Alleles Revealed 653 by Controlled Sexual Recombination in Pea Aphids. Genetics. 2015;199:581–93. 654 28. Stewart SA, Hodge S, Ismail N, Mansfield JW, Feys BJ, Prospéri J-M, et al. The RAP1 Gene Confers 655 Effective, Race-Specific Resistance to the Pea Aphid in Medicago truncatula Independent of the 656 Hypersensitive Reaction. Mol Plant-Microbe Interactions®. 2009;22:1645–55. 657 29. Julier B, Huguet T, Chardon F, Ayadi R, Pierre J-B, Prosperi J-M, et al. Identification of quantitative 658 trait loci influencing aerial morphogenesis in the model legume Medicago truncatula. Theor Appl 659 Genet. 2007;114:1391–406. 660 30. Peccoud J, Mahéo F, de la Huerta M, Laurence C, Simon J-C. Genetic characterisation of new 661 host-specialised biotypes and novel associations with bacterial symbionts in the pea aphid complex. 662 Insect Conserv Divers. 2015;8:484–92. 663 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 21 31. Carolan JC, Caragea D, Reardon KT, Mutti NS, Dittmer N, Pappan K, et al. Predicted Effector 664 Molecules in the Salivary Secretome of the Pea Aphid ( Acyrthosiphon pisum ): A Dual 665 Transcriptomic/Proteomic Approach. J Proteome Res. 2011;10:1505–18. 666 32. Carolan JC, Fitzroy CIJ, Ashton PD, Douglas AE, Wilkinson TL. The secreted salivary proteome of 667 the pea aphid Acyrthosiphon pisum characterised by mass spectrometry. PROTEOMICS. 668 2009;9:2457–67. 669 33. Pitino M, Hogenhout SA. Aphid Protein Effectors Promote Aphid Colonization in a Plant Species-670 Specific Manner. Mol Plant-Microbe Interactions®. 2013;26:130–9. 671 34. Kaloshian I, Walling LL. Hemipterans as Plant Pathogens. Annu Rev Phytopathol. 2005;43:491–672 521. 673 35. Furch ACU, Van Bel AJE, Will T. Aphid salivary proteases are capable of degrading sieve-tube 674 proteins. J Exp Bot. 2015;66:533–9. 675 36. Cantón PE, Bonning BC. Extraoral digestion: outsourcing the role of the hemipteran midgut. Curr 676 Opin Insect Sci. 2020;41:86–91. 677 37. Ferrari J, Via S, Godfray HCJ. POPULATION DIFFERENTIATION AND GENETIC VARIATION IN 678 PERFORMANCE ON EIGHT HOSTS IN THE PEA APHID COMPLEX. Evolution. 2008;62:2508–24. 679 38. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-680 seq aligner. Bioinformatics. 2013;29:15–21. 681 39. Thorpe P, Escudero-Martinez CM, Cock PJA, Eves-van Den Akker S, Bos JIB. Shared 682 Transcriptional Control and Disparate Gain and Loss of Aphid Parasitism Genes. Genome Biol Evol. 683 2018;10:2716–33. 684 40. Thorpe P. Genome annotations for: Multi-omics approaches define novel aphid effector 685 candidates associated with virulence and avirulence phenotypes. 2024. 686 41. Prosser WA, Douglas AE. The aposymbiotic aphid: An analysis of chlortetracycline-treated pea 687 aphid, Acyrthosiphon pisum. J Insect Physiol. 1991;37:713–9. 688 42. Petersen TN, Brunak S, von Heijne G, Nielsen H. SignalP 4.0: discriminating signal peptides from 689 transmembrane regions. Nat Methods. 2011;8:785–6. 690 43. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting 691 for selection bias. Genome Biol. 2010;11:R14. 692 44. Perez-Riverol Y, Bai J, Bandla C, García-Seisdedos D, Hewapathirana S, Kamatchinathan S, et al. 693 The PRIDE database resources in 2022: a hub for mass spectrometry-based proteomics evidences. 694 Nucleic Acids Res. 2022;50:D543–52. 695 696 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 22 Tables 697 Table 1 Genes and proteins overlapping in multiple experiments. All genes shown that are 698 represented in at least three datasets, plus all genes intersected between F1 transcriptome and at 699 least one other dataset. Saliva and salivary gland data are proteins, head and body data are 700 transcripts. A. Proteins and upregulated genes in virulent aphids (N116, VIR F1 pool); B. Proteins and 701 upregulated genes in avirulent aphids (PS01, AVR F1 pool). Y = protein present and/or RNA 702 differentially expressed. Full gene and protein lists are in Supplementary Material 3 and 4. 703 Gene Annotation Saliva Salivary gland Parent head Parent body F1 body A) N116 & VIR F1 ACPISUM_000319 ACYPI007553 Y Y Y ACPISUM_006458 aldo-keto reductase family 1 member B10-like Y Y Y ACPISUM_025240 aminopeptidase N Y Y Y Y ACPISUM_005699 aminopeptidase N Y Y Y Y ACPISUM_025168 aminopeptidase N Y Y Y Y ACPISUM_009258 aminopeptidase N Y Y Y ACPISUM_024778 aminopeptidase N Y Y Y ACPISUM_026844 aminopeptidase N Y Y Y ACPISUM_025015 aminopeptidase N Y Y Y Y ACPISUM_023906 Apoptosis inducing protein Y Y Y ACPISUM_020864 F-actin-capping protein subunit alpha Y Y Y ACPISUM_023535 glutamate-gated chloride channel-like Y Y ACPISUM_019971 glutathione hydrolase 1 proenzyme-like Y Y Y ACPISUM_010531 hypothetical protein X975_16721 Y Y ACPISUM_013751 LYR motif-containing protein 4 Y Y Y ACPISUM_013796 myrosinase 1-like Y Y ACPISUM_006164 ---NA--- Y Y Y ACPISUM_023321 papain inhibitor-like Y Y Y ACPISUM_009624 proline-rich extensin-like protein EPR1 Y Y Y ACPISUM_028519 single-stranded DNA-binding replication protein A Y Y Y ACPISUM_025560 ubiquinone biosynthesis monooxygenase COQ6, mitochondrial Y Y Y ACPISUM_008675 uncharacterized protein LOC100162547 Y Y Y ACPISUM_007320 uncharacterized protein LOC100167449 Y Y Y Y ACPISUM_001031 uncharacterized protein LOC100571631 Y Y Y ACPISUM_016519 uncharacterized protein LOC100573156 Y Y Y ACPISUM_010687 uncharacterized protein LOC103309122 Y Y Y ACPISUM_017388 uncharacterized protein LOC103309964 Y Y Y Y ACPISUM_009099 uncharacterized protein LOC112598674 Y Y Y ACPISUM_027918 vacuolar protein sorting-associated protein 29 Y Y Y B) PS01 & AVR F1 ACPISUM_000957 AGAP002382-PA-like protein Y Y Y ACPISUM_015173 AGAP011571-PA-like protein Y Y ACPISUM_002223 aminopeptidase N Y Y Y ACPISUM_003737 aminopeptidase N Y Y Y Y ACPISUM_023448 aminopeptidase N Y Y Y Y ACPISUM_028967 aminopeptidase N Y Y Y Y ACPISUM_021545 aminopeptidase N Y Y Y ACPISUM_009259 aminopeptidase N Y Y Y ACPISUM_009580 anoctamin-1-like Y Y Y ACPISUM_012705 CD63 antigen Y Y Y ACPISUM_006933 cuticular protein Y Y Y ACPISUM_019160 glutathione S-transferase 1-1-like Y Y Y Y ACPISUM_019168 glutathione S-transferase 1-1-like Y Y Y Y ACPISUM_001883 glutathione S-transferase D7-like Y Y Y ACPISUM_016389 histone acetyltransferase KAT6B isoform X1 Y Y ACPISUM_009097 multidrug resistance-associated protein 1 Y Y Y ACPISUM_011553 ---NA--- Y Y Y ACPISUM_011754 ---NA--- Y Y Y Y ACPISUM_004702 ---NA--- Y Y ACPISUM_021569 ---NA--- Y Y .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 23 ACPISUM_025236 ---NA--- Y Y ACPISUM_014327 ---NA--- Y Y Y ACPISUM_016390 ---NA--- Y Y Y ACPISUM_017200 ---NA--- Y Y Y ACPISUM_027631 ---NA--- Y Y Y ACPISUM_028853 ---NA--- Y Y Y ACPISUM_019381 neural cell adhesion molecule L1 isoform X1 Y Y Y ACPISUM_020816 peroxidase-like Y Y Y Y ACPISUM_019857 peroxidase-like Y Y Y ACPISUM_019870 peroxidasin homolog Y Y Y ACPISUM_000958 phospholipase DDHD2-like Y Y ACPISUM_006758 piggyBac transposable element-derived protein 4-like Y Y Y ACPISUM_022113 piwi-like protein Siwi Y Y ACPISUM_010778 predicted protein Y Y Y ACPISUM_019013 protein ABHD18 Y Y Y ACPISUM_021997 regucalcin-like Y Y Y Y ACPISUM_021999 regucalcin-like Y Y Y ACPISUM_001383 replication protein A 70 kDa DNA-binding subunit-like Y Y Y ACPISUM_015166 TBC1 domain family member 19 Y Y Y ACPISUM_014232 tubulin glycylase 3A-like Y Y ACPISUM_008377 uncharacterized family 31 glucosidase KIAA1161-like Y Y Y ACPISUM_008379 uncharacterized family 31 glucosidase KIAA1161-like Y Y Y ACPISUM_008380 uncharacterized family 31 glucosidase KIAA1161-like Y Y Y ACPISUM_012348 uncharacterized protein LOC100158692 Y Y Y ACPISUM_018433 uncharacterized protein LOC100158721 Y Y ACPISUM_007487 uncharacterized protein LOC100160601 Y Y Y ACPISUM_016065 uncharacterized protein LOC100161530 Y Y Y ACPISUM_007076 uncharacterized protein LOC100163035 Y Y Y ACPISUM_016064 uncharacterized protein LOC100570074 Y Y Y ACPISUM_029311 uncharacterized protein LOC100570454 Y Y Y Y ACPISUM_008664 uncharacterized protein LOC100570454 Y Y Y ACPISUM_007394 uncharacterized protein LOC100572241 Y Y Y Y ACPISUM_021703 uncharacterized protein LOC100575642 Y Y Y ACPISUM_029930 uncharacterized protein LOC100575698 Y Y Y Y Y ACPISUM_006906 uncharacterized protein LOC100575848 Y Y Y ACPISUM_003989 uncharacterized protein LOC103307823 Y Y Y ACPISUM_024374 uncharacterized protein LOC107882950 Y Y ACPISUM_015285 uncharacterized protein LOC107883982 Y Y Y ACPISUM_000491 uncharacterized protein LOC111028731 Y Y Y ACPISUM_027814 uncharacterized SDCCAG3 family protein-like Y Y Y 704 705 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 24 Table 2 Significantly enriched GO terms within differentially expressed transcript and 706 protein data. Terms enriched at FDR<0.05, after manual curation to remove redundancies, 707 retaining the terms with lowest FDR. Full lists of enriched terms are in Supplementary 708

5. 709 GO Category Term Ontology group No. genes in DE set P value FDR Whole body RNA N116 up 0004177 aminopeptidase activity MF 8 2.72E-06 0.0270 PS01 up 0004177 aminopeptidase activity MF 9 1.43E-07 0.0014 0017177 glucosidase II complex CC 4 2.09E-06 0.0104 AVR F1 up 0017177 glucosidase II complex CC 3 2.42E-06 0.0240 VIR F1 up No terms Heads RNA N116 up on A17 0045271 respiratory chain complex I CC 11 1.55E-07 9.63E-05 0005743 mitochondrial inner membrane CC 22 2.01E-07 0.0001 0004177 aminopeptidase activity MF 10 2.13E-06 0.0009 0016491 oxidoreductase activity MF 31 3.24E-06 0.0012 0005875 microtubule associated complex CC 25 1.72E-05 0.0052 0019395 fatty acid oxidation BP 5 2.44E-05 0.0064 0042826 histone deacetylase binding MF 4 0.00011 0.0275 0045239 tricarboxylic acid cycle enzyme complex CC 3 0.00014 0.0327 0004448 isocitrate dehydrogenase activity MF 3 0.00015 0.0338 N116 up on DZA 0006635 fatty acid beta-oxidation BP 5 1.05E-06 0.0082 0004177 aminopeptidase activity MF 9 2.66E-06 0.0082 0004449 isocitrate dehydrogenase (NAD+) activity MF 3 1.40E-05 0.0198 0006099 tricarboxylic acid cycle BP 6 3.54E-05 0.0389 PS01 up on A17 0004177 aminopeptidase activity MF 10 7.77E-07 0.00771 PS01 up on DZA 0004177 aminopeptidase activity MF 11 5.59E-09 5.53E-05 Salivary gland proteins N116 up 0003983 UTP:glucose-1-phosphate uridylyltransferase activity MF 3 4.97E-06 0.0213 PS01 up No terms 710 711 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 25 Table 3. Comparison of expression patterns of exopeptidases detected in saliva and salivary 712 glands. All detected proteins are listed, along with whether they were differentially expressed, and 713 whether the patterns were also reflected in the transcriptomes. Sal = saliva; SG = salivary gland; Y = 714 protein present. 715 716 717 ACPISUM v3 gene Nearest ACYPI gene(s) Sal Sal SG SG Sal SG heads whole N116 PS01 N116 PS01 Aminopeptidase N ACPISUM_005699-T1 ACYPI080623 ACYPI070600 ACYPI005810 Y Y Y Y ACPISUM_025168-T1 ACYPI068031 Y Y ACPISUM_027632-T1 ACYPI073645 Y Y ACPISUM_009259-T1 ACYPI007868 Y Y Y ACPISUM_009258-T1 ACYPI007868 Y Y ACPISUM_024778-T1 ACYPI072916 Y ACPISUM_025015-T1 ACYPI061522 ACYPI21510 Y Y ACPISUM_023174-T1 ACYPI006366 Y Y ACPISUM_025240-T1 ACYPI070333 Y Y Y ACPISUM_000115-T1 ACYPI000001 Y ACPISUM_003737-T1 ACYPI067691 Y Y Y ACPISUM_021545-T1 ACYPI085147 ACYPI002583 Y Y Y Y ACPISUM_029674-T1 ACYPI072988 Y Y Y Y ACPISUM_023448-T1 ACYPI010198 Y Y Y Y ACPISUM_028967-T1 ACYPI083965 Y Y Y Y ACPISUM_000246-T1 ACYPI086097 ACYPI43770 ACYPI068046 Y Y ACPISUM_010796-T1 ACYPI071232 ACYPI33244 Y Y Y Y ACPISUM_012062-T1 ACYPI22813 Y Y Y Y ACPISUM_018507-T1 ACYPI21711 ACYPI084528 ACYPI003165 Y Y Y Y ACPISUM_006298-T1 ACYPI41708 ACYPI22605 Y Y Y Y ACPISUM_019635-T1 ACYPI060722 Y ACPISUM_019937-T1 ACYPI49161 Y Y Y Y ACPISUM_026119-T1 ACYPI083984 Y Y Y Y ACPISUM_017858-T1 ACYPI54528 ACYPI001911 Y Y ACPISUM_002219-T1 ACYPI44040 Y Y Y Y ACPISUM_014203-T1 ACYPI067721 Y Y ACPISUM_018506-T1 ACYPI21557 Y Y ACPISUM_019609-T1 ACYPI001203 Y Y ACPISUM_019610-T1 ACYPI071951 Y Y Aminopeptidase N total proteins detected 24 20 21 17 Angiotensin converting enzyme ACPISUM_008374-T1 ACYPI000733 Y Y Y Y ACPISUM_024301-T1 ACYPI084554 Y Y Y Y ACPISUM_024303-T1 ACYPI071320 Y Y Y Y ACPISUM_020790-T1 ACYPI008911 Y Y Key not detected detected, not DE up in N116 up in PS01 Presence/absence Differential expression .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 26 Figure legends 718 Figure 1. Summary of transcriptome and proteome analysis pipeline. Resistant and susceptible host 719 plants carried or lacked the RAP1 aphid resistance QTL, respectively. For all experiments, virulent N116 720 and avirulent PS01 aphids were compared. In addition, BSA-RNA-Seq was done on whole body pooled 721 samples of F1 virulent and avirulent aphids. 722 Figure 2. Virulence phenotypes of parental aphid clones and selections from F1 populations used 723 for BSA-RNA-Seq. Tested on two M. truncatula genotypes carrying the RAP1 locus: Jemalong A17 and 724 a resistant near isogenic line (RNIL) derived from a cross between A17 and DZA315.16. The parental 725 genotypes and selections from the F1 populations shown here were all used for the BSA -RNA-Seq 726 experiment. Data are expressed as virulence index, assessed 10 d after infestation. Phenotypes of F1 727 clones were classified using the following virulence index cut -offs: A17 VIR >4, AVR 4, 728 AVR <4.5. Orange circles are NP (N116 female x PS01 male); blue triangles are PN (PS01 female x N116 729 male); red is N116, and green is PS01, with each of three parental data points from a separate batch 730 of F1 tests. The full population phenotype data are provided in Supplementary Material 1. 731 Figure 3. Transcriptome analysis of aphid heads. Samples were dissected heads from PS01 and N116 732 genotypes infested on Medicago truncatula A17 or DZA315.16 for 24 h, with n=3 biological replicates. 733 Aphid genotype PS01 is avirulent on M. truncatula A17, and all other combinations represent 734 compatible interactions. A. Clustering of transcriptional responses of pea aphid, showing samples 735 clustered more strongly based on aphid genotype than on host interaction; B. Principal components 736 analysis. The top two principal components explain >68% of the variation among transcriptional 737 responses. Samples group largely by aphid genotype rather than host interaction; C. Numbers of genes 738 differentially expressed between the different aphid genotypes on different M. truncatula genotypes. 739 Of the 935 DE genes between PS01 and N116 on A17, 483 were up in N116 and 452 were up in PS01. 740 Of the 758 DE genes on DZA hosts, 395 were up in N116 and 363 were up in PS01. Accompanying gene 741 lists and annotations are provided in Supplementary Material 3. 742 Figure 4. Transcriptome analysis of whole aphids. Aphids were infested on Medicago truncatula A17 743 for 24 h, with n=5 biological replicates. A . Clustering of transcriptional responses of pea aphid PS01, 744 N116, bulked F1 VIR and AVR progeny. Responses within biological replicates are more strongly 745 correlated than responses among different aphid genotypes; B. Principal components analysis. The 746 top 2 principal components explain >45% of the variation among transcriptional responses of PS01, 747 N116, AVR and VIR F1 progeny replicates, and separate the responses of the different aphid genotypes 748 and F1 pools. C. Numbers of genes differentially expressed between the different aphid genotypes 749 and pools. Accompanying gene lists and annotations are provided in Supplementary Material 3. 750 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 27 Figure 5. Differential gene expression in pea aphid genotypes N116, PS01 and bulked F1 pools of 751 virulent and avirulent progeny. A. Numbers of genes up - versus down-regulated in comparisons of 752 parent genotypes N116 (virulent) and PS01 (avirulent), and their F1 progeny pools, all on A17 host 753 plants. Orange bars represent the numbers of genes up-regulated in genotype N116 or the VIR F1 pool 754 compared with the genotype PS01 and the AVR pool, respectively. Blue bars represent the numbers 755 of genes up-regulated in genotype PS01 or the AVR F1 pool compared with the genotype N116 and 756 the VIR pool, respectively; B. Overlaps in genes up-regulated in genotype N116 whole body and head 757 tissues compared to genotype PS01, and up-regulated in the VIR F1 pool compared to the AVR F1 pool; 758 C. Overlaps in genes up -regulated in genotype PS01 whole body and head tissues compared to 759 genotype N116, and up-regulated in the AVR F1 pool compared to the VIR F1 pool; D. Overlaps in up-760 regulated genes among whole body transcriptomes of N116, PS01, VIR F1 pool and AVR F1 pool. 761 Figure 6. Comparative proteomic analysis of salivary glands and saliva for pea aphid genotypes N116 762 and PS01. Venn diagrams of the number of proteins shared and found exclusively in A) salivary glands 763 and B) saliva identified for both genotypes. Principal Components Analysis (PCA) of C) salivary glands 764 and D) saliva distinguishes both genotypes clearly. Volcano plots based on -log10 p values and log2 fold 765 differences highlighting the statistically significant differentially abundant (SSDA) proteins (p≤0.05) for 766 E) salivary glands and F) saliva. Annotations are shown for the top 12 proteins of increased and 767 decreased abundances. 768 Figure 7. Selected differentially expressed genes from whole body transcriptomes . A,B 769 representative genes upregulated both in virulent parent and in virulent F1 pool; C -F representative 770 genes upregulated both in avirulent parent and in avirulent F1 pool. G,H representative genes with 771 opposite regulation between parent and F1 pairs. Each point represents an individual RNA-Seq library 772 (n=5). *** indicates FDR<0.001. 773 774 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 28

775 Funding 776 We thank the Biotechnology and Biological Sciences Research Council for funding to CT 777 (BB/N002830/1) and JB (BB/N002660/1). We thank Umer Rashid and Martin Selby for expert 778 technical assistance. 779 Author contributions 780 JB, CT, JC, PT, SA and RLC designed the experiments. PT, SA, RLC, ND, JCS, JI and SK conducted the 781 experiments. PT, SA, RLC, JB, CT, JC, ND and JCS analysed the data. CT, JB, JC and PT wrote the paper. 782 All authors approved the submitted manuscript. 783 Conflicts 784 The authors declare that they have no competing interests. 785 Supplementary Materials 786 Supplementary Material 1. F1 aphid phenotyping. 787 Supplementary Material 2. Read mapping summary. 788 Supplementary Material 3. RNA-Seq Head and whole body differentially expressed genes. 789 Supplementary Material 4. Salivary gland and saliva proteomics. 790 Supplementary Material 5. GO enrichment791 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 29 792 793 Figure 1. Summary of transcriptome and proteome analysis pipeline. Resistant and susceptible host plants carried or lacked the RAP1 aphid resistance QTL, 794 respectively. For all experiments, virulent N116 and avirulent PS01 aphids were compared. In addition, BSA-RNA-Seq was done on whole body pooled samples 795 of F1 virulent and avirulent aphids. 796 797 Virulent vs avirulent aphids Differential protein abundance Proteomics Heads RNA-Seq Salivary gland SalivaWhole body Differential gene expression Resistant vs susceptible host Candidate Effectors .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 30 798 799 Figure 2. Virulence phenotypes of parental aphid clones and selections from F1 populations used for BSA-RNA-Seq. Tested on two M. truncatula genotypes 800 carrying 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 and 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 infestation. Phenotypes of F1 clones were classified using the following virulence index cut -offs: A17 VIR >4, AVR 4, AVR <4.5. Orange circles 803 are 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 from a separate batch of F1 tests. The full population phenotype data are provided in Supplementary Material 1. 805 806 0 1 2 3 4 5 6 8 9 0 1 2 3 4 5 6 8 9 Virulence index on rNIL Virulence index on A1 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 31 807 Figure 3. Transcriptome analysis of aphid heads. Samples were dissected heads from PS01 and N116 genotypes infested on Medicago truncatula A17 or 808 DZA315.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 interactions. A. Clustering of transcriptional responses of pea aphid, showing samples clustered more strongly based on aphid genotype than on host 810 interaction; B. Principal components analysis. The top two principal components explain >68% of the variation among transcrip tional responses. Samples 811 group largely by aphid genotype rather than host interaction; C. Numbers of genes differentially expressed between the different aphid genotypes on different 812 M. 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 hosts, 395 were up in N116 and 363 were up in PS01. Accompanying gene lists and annotations are provided in Supplementary Material 3. 814 A B PS01/DZAPS01/A17N116/DZAN116/A17 845935330N116/A17 758782033N116/DZA 290782935PS01/A17 029758845PS01/DZA C .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 32 815 Figure 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 transcriptional responses of pea aphid PS01, N116, bulked F1 VIR and AVR progeny. Responses within biological replicates are more strongly correlated than 817 responses among different aphid genotypes; B. Principal components analysis. The top 2 principal components explain >45% of t he variation among 818 transcriptional responses of PS01, N116, AVR and VIR F1 progeny replicates, and separate the responses of the different aphid genotypes and F1 pools. C. 819 Numbers of genes differentially expressed between the different aphid genotypes and pools. Accompanying gene lists and annota tions are provided in 820 Supplementary Material 3. 821 822 B Avr F1 poolPS01Vir F1 poolN116 2695771480N116 882410148Vir F1 pool 2460241577PS01 024688269Avr F1 pool C A .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 33 823 Figure 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 versus 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 bars 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 Blue 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 B. 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 to 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 AVR 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 132 165 317 21 2 1 N116 up - heads N116 up - whole aphid Vir F1 pool - up B 260 168 79 31 9 22 2 PS01 up - whole aphid PS01 up - heads Avr F1 pool - up C -500 -400 -300 -200 -100 0 100 200 300 400 500 A 299 24 452 64 number of DE genes avirulentvirulent 278 483 D N116 up - whole aphid PS01 up - whole aphid Avr F1 pool - up Vir F1 pool -up 294 234 31 139 30 2 3 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 34 831 Figure 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 both 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 differentially 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 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint 35 835 Figure 7. Selected differentially expressed genes from whole body transcriptomes . A,B representative genes upregulated both in virulent parent and in 836 virulent F1 pool; C-F representative genes upregulated both in avirulent parent and in avirulent F1 pool. G,H representative genes with opposite regulation 837 between parent and F1 pairs. Each point represents an individual RNA-Seq library (n=5). *** indicates FDR<0.001. 838 ‐2 ‐1 0 1 2 3 ero centred log2 expression N116 PS01 F1 VIR F1 AVR ACPIS M 006933 cu cular protein ‐3 ‐2 ‐1 0 1 2 3 ero centred log2 expression N116 PS01 F1 VIR F1 AVR ACPIS M 0198 0 peroxidasin homolog ‐2 ‐1 0 1 2 ero centred log2 expression N116 PS01 F1 VIR F1 AVR ACPIS M 005464 transcrip onal regulator ERG homolog ‐3 ‐2 ‐1 0 1 2 3 ero centred log2 expression N116 PS01 F1 VIR F1 AVR ACPIS M 013 96 myrosinase 1‐like ‐2 ‐1 0 1 2 ero centred log2 expression N116 PS01 F1 VIR F1 AVR ACPIS M 0199 1 glutathione hydrolase 1 proenzyme‐like ‐4 ‐3 ‐2 ‐1 0 1 2 3 ero centred log2 expression N116 PS01 F1 VIR F1 AVR ACPIS M 008380 family 31 glucosidase KIAA1161‐like ‐2 ‐1 0 1 2 ero centred log2 expression N116 PS01 F1 VIR F1 AVR ACPIS M Api15 8 uncharacterisedprotein ‐3 ‐2 ‐1 0 1 2 3 ero centred log2 expression N116 PS01 F1 VIR F1 AVR ACPIS M 02199 regucalcin‐like A B C D E F G H ACPISUM_029930 .CC-BY 4.0 International licensemade available under a (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 The copyright holder for this preprintthis version posted July 30, 2024. ; https://doi.org/10.1101/2024.07.30.605808doi: bioRxiv preprint