Strong genotype‐by‐genotype interactions between aphid‐defensive symbionts and parasitoids persist across different biotic environments

Abstract The dynamics of coevolution between hosts and parasites are influenced by their genetic interactions. Highly specific interactions, where the outcome of an infection depends on the precise combination of host and parasite genotypes (G × G interactions), have the potential to maintain genetic variation by inducing negative frequency‐dependent selection. The importance of this effect also rests on whether such interactions are consistent across different environments or modified by environmental variation (G × G × E interaction). In the black bean aphid, Aphis fabae, resistance to its parasitoid Lysiphlebus fabarum is largely determined by the possession of a heritable bacterial endosymbiont, Hamiltonella defensa, with strong G × G interactions between H. defensa and L. fabarum. A key environmental factor in this system is the host plant on which the aphid feeds. Here, we exposed genetically identical aphids harbouring three different strains of H. defensa to three asexual genotypes of L. fabarum and measured parasitism success on three common host plants of A. fabae, namely Vicia faba, Chenopodium album and Beta vulgaris. As expected, we observed the pervasive G × G interaction between H. defensa and L. fabarum, but despite strong main effects of the host plants on average rates of parasitism, this interaction was not altered significantly by the host plant environment (no G × G × E interaction). The symbiont‐conferred specificity of resistance is thus likely to mediate the coevolution of A. fabae and L. fabarum, even when played out across diverse host plants of the aphid.

to counter-adapt and is soon going to be at a disadvantage compared to rare genotypes. Many rare genotypes are therefore favoured over few frequent ones in the long term (Clarke, 1976;Judson, 1995). G × G interactions may also contribute to the evolutionary maintenance of sexual reproduction and recombination, which promote diversity by combining existing genes and producing new, rare offspring genotypes (The Red Queen Hypothesis: Bell, 1982;Hamilton, 1980;Hamilton et al., 1990;Jaenike, 1978).
The potentially far-reaching influence of G × G interactions on host-parasite coevolution also hinges on their sensitivity to environmental variation. With the environment varying across the geographic distribution of interacting species, also the conditions for successful host infection or defence can change, potentially altering the strength and direction of selection acting on host and parasite (Nuismer et al., 2000;Tétard-Jones et al., 2007;Thompson, 2005;Wolinska & King, 2009). Different genotypes of a species may respond unequally to changing environmental conditions, resulting in genotype-by-environment (G × E) interactions that influence the outcome of host-parasite encounters. Moreover, the strength and specificity of G × G interactions may be dependent on the environment, which is referred to as a genotype-by-genotype-by-environment interaction (G × G × E). Such three-way interactions may contribute to the maintenance of genetic diversity at the species level as well, and they make the outcome of genotype-specific selection of hosts and parasites less predictable (Mostowy & Engelstädter, 2011;Wolinska & King, 2009). Experimental evidence for G × G interactions that are sensitive to environmental variation, that is for G × G × E interactions, is not abundant but exists for diverse systems. Examples include G × G × E interactions with the rhizosphere environment influencing the outcome of specific interactions between plants and herbivore genotypes (Tétard-Jones et al., 2007), or the food environment influencing the G × G interactions that determine parasite infection in bumblebees (Sadd, 2011). Similarly, Bryner and Rigling (2011) found that temperature interacts with the genotypes of tree pathogenic fungi and their hyperparasitic viruses when predicting fungal virulence, and a three-way interaction between temperature, host genotype and parasite genotype may determine the transmission potential of viral diseases by mosquito vectors (Zouache et al., 2014). Apart from experimental evidence, modelling approaches support the potential influence of environmental changes on hostparasite coevolution through three-way interactions (Mostowy & Engelstädter, 2011). Taken together, these references underline the importance of incorporating environmental variability into classical G × G interaction studies, prior to generalizing conclusions to more complex natural systems.
Aphids have become popular model organisms for studying host-parasite interactions, not least because of the fascinating way by which certain species resist natural enemies: they carry defensive symbionts protecting them, for example, against pathogenic fungi (Lukasik et al., 2013;Scarborough et al., 2005) or parasitoid wasps Oliver et al., 2003;Vorburger et al., 2009). The best-studied example for the latter is the gammaproteobacterium Hamiltonella defensa (Moran, Russell, et al., 2005;Oliver et al., 2003), a maternally transmitted endosymbiont. There are different strains of the endosymbiont H. defensa which may confer different levels of resistance against parasitoid wasps, depending on the species or also the genotype of the attacking parasitoids Cayetano & Vorburger, 2015). Because maternal transmission of H. defensa is very reliable (Darby & Douglas, 2003;Peccoud et al., 2014;Vorburger et al., 2017), H. defensa may be regarded as a form of a selectable resistance trait of its aphid host, with different strains acting as different genotypes determining the characteristics of this trait (Jaenike, 2012). Since H. defensa-conferred resistance tends to be much stronger than the basal resistance of the aphid (Oliver et al., 2005;Vorburger et al., 2009), coevolution between aphid hosts and parasitoids is likely mediated by defensive symbionts (Vorburger & Perlman, 2018).
The present study is concerned with the interaction between the black bean aphid, Aphis fabae (Hemiptera: Aphididae), and its main parasitoid Lysiphlebus fabarum (Hymenoptera: Braconidae, Aphidiinae). Different experiments have shown that the rate of successful parasitism in this system is determined by strong G × G interactions between L. fabarum and the aphids' symbiont H. defensa (Rouchet & Vorburger, 2012;Schmid et al., 2012). This suggests a decisive role of H. defensa in governing natural coevolutionary dynamics between A. fabae and L. fabarum (Kwiatkowski et al., 2012).
However, many of the experiments yielding the current knowledge were justifiably done under constant ambient conditions in the laboratory, which leaves open the question of how consistent-and thus relevant-such G × G interactions may be in heterogeneous natural environments. The only environmental variable that has been explicitly manipulated in this system is temperature: on the one hand, Cayetano and Vorburger (2013b) showed that G × G interactions between L. fabarum and H. defensa remain qualitatively the same at different ambient temperatures, even though the level of resistance conferred by H. defensa drops with increasing temperature, as also seen in pea aphids (Bensadia et al., 2006;Doremus et al., 2018). On the other hand, the same authors found that heat shocks experienced by the aphids could affect G × G interactions, albeit only at very high temperatures that are rarely experienced in nature (Cayetano & Vorburger, 2013a). We aimed to complement these studies by manipulating a different, yet crucially important environmental variable of the same host-endosymbiont-parasitoid system, namely the aphid's host plant.  McLean et al. (2011) to suggest that selection pressure by parasitoids varies across different host plants. As shown for the same aphid species, also predation rates can be influenced by the host plant (Aquilino et al., 2005). Generally, host plant choice may affect the fitness of herbivorous insects and consequently the fitness of their predators or parasitoids (Pan et al., 2020).
While the influence of the host plant on insect interactions is striking in many systems, actual studies of G × G × E interactions with host plant as the environmental variable are rare. Investigating such interactions requires a study system where specific genotype combinations can be replicated, as is the case for the A. fabae/H. defensa/L. fabarum system. A. fabae can reproduce clonally, and its infection with H. defensa can readily be manipulated by microinjection (Oliver et al., 2003;Sochard et al., 2020), while the occurrence of asexual reproduction in L. fabarum (thelytoky, see Sandrock et al. (2011)) enables the use of distinct, genetically homogeneous lines also for the parasitoid. Environmental variability due to host plants is of high relevance for the multivoltine Aphis fabae, where opportunistic switches between crops and weeds from one aphid generation to the other are important to allow continuous feeding and reproduction over a whole growing season. We thus investigated the influence of three different, common host plants of A. fabae on G × G interactions between the aphid-defensive symbiont H. defensa and the parasitoid L. fabarum, using a full factorial design. The host plant had a strong effect on overall parasitism rates and thus on wasp reproductive success, and it also affected aphid fitness independent of parasitism rates. The host plant did, however, not alter the strong G × G specificity between L. fabarum and H. defensa. Hence, our results support earlier studies, suggesting that in our model system, coevolution between aphids and parasitoids is largely symbiontmediated and governed by genotype-specific interactions, which remain remarkably stable across different environments.
We chose four aphid lines from our laboratory collection showing a broad spectrum of resistance to different lines of the parasitoid L. fabarum, as seen in previous studies (Cayetano & Vorburger, 2013a, 2013b, 2015. The four lines originate from a single clone of A. fabae fabae (clone ID: 407), which was collected in Switzerland in 2006 from Chenopodium. One of these lines was free of any facultative symbionts (407), and three lines were uniquely infected with one of three genetically different strains of H. defensa: H15, H76 and H402. To obtain these lines, the H. defensa-free aphid clone 407 had been infected by microinjection of haemolymph from H. defensacarrying aphid clones (Cayetano & Vorburger, 2015), resulting in stable, heritable infections. The infected lines are referred to as 407 H15 , 407 H76 and 407 H402 and have been maintained parthenogenetically since the infection. Their identity was reconfirmed immediately before the experiment by microsatellite genotyping of the aphid (Coeur d'acier et al., 2004) and sequencing of the H. defensa gene murE (Degnan & Moran, 2008b). Using a single aphid clone should exclude genetic variation beyond the endosymbiont strain in the aphid lines.
We included the H. defensa-free aphid line in order to relate levels of H. defensa-conferred resistance to the aphid's basal resistance.
As parasitoids, we used three asexual, isofemale lines of L. fabarum (06-242, 07-64 and 09-369), the most frequent parasitoid species of A. fabae in Switzerland (Rothacher et al., 2016). The lines had been started from single asexual females collected between 2006 and 2009. Parasitoid wasps oviposit single eggs into aphids. After hatching, the wasp larva feeds on the aphid's body, eventually kills the aphid and pupates within the emptied aphid exoskeleton. At this stage, a parasitized aphid is clearly recognizable and referred to as a mummy. replicate helped to level out potential environmental effects carried over from the stock cultures and allowed for physiological adaptation to the different plant species. To initiate the experimental generation, we put four adult female aphids on a plant on day 1. We let them reproduce for up to two days in order to have similar numbers of nymphs on each plant, before removing the adults and counting the nymphs on day 3. On day 4, we added two female parasitoids to each plant and removed them again after six hours. We then waited until the aphid mummies (= parasitized aphids) were clearly visible on day 13 and counted the parasitized as well as the adult non-parasitized aphids. As response variable, we used the number of mummies divided by the number of nymphs initially exposed to the wasps (parasitism rate). A considerable number of aphids died in the time between the counting of the aphid nymphs and the counting of the mummies. To ensure that this would not falsify our conclusions, we did a second analysis with parasitism rate calculated as the number of mummies divided by the number of aphids still present (alive or parasitized) on the plant when counting mummies.

| Experimental set-up
We further measured the fresh weight of the mothers of the experimental aphid generation on a precision balance (MX5, Mettler Toledo, Greifensee, Switzerland) to assess potential effects of host plant and aphid line on aphid size. As a measure connected to parasitoid fitness, we determined the proportion of wasps hatching from the mummies we collected during the experiment (emergence rate).
For this, we cut the whole plant after having counted the mummies and kept it in an air permeable bag until the adult wasps emerged from the mummies.
We performed a three-way factorial analysis of deviance testing for the effects of aphid line, parasitoid line and host plant as well as the two-and three-way interactions, and for the main effect of experimental block. We used the function Anova from the R package car v3.0.7 (Fox & Weisberg, 2019) with F tests as recommended by Crawley (2014) for quasilikelihood fits. We treated experimental block as fixed since quasibinomial errors are not implemented for generalized mixed models in R. For consistency, we treated block as a fixed effect in all further analyses. While the initial number of nymphs exposed to the wasps differed between plants (Beta 15.6 ± 7.2 SD, Chenopodium 13.8 ± 6.0, Vicia: 17.3 ± 6.7), we did not include this value in the final model, since it had no significant effect on parasitism rates when aphid line and host plant were included (Table S1), suggesting that parasitoids were host limited (and not time limited) in our assays. We did the analysis once for the full data set (all aphid lines) and once with the data set restricted to the H. defensa-infected aphid lines. Only in the latter case does the aphid line × parasitoid line interaction strictly reflect the G × G interactions between symbionts and parasitoids. Certain treatment combinations resulted in zero mummies in all replicates, and thus a group variance of zero, which led to problems with model convergence. To avoid this, we edited our data such that we manually added one mummy to one replicate of each 'zero parasitism' treatment combination. This minor intervention should have reduced treatment differences and hence made comparisons more conservative.
The parasitoid emergence rate was analysed with generalized linear models with logit link functions and binomial errors. Since we could analyse only samples where at least one mummy had formed, we performed separate analyses on data subsets from the aphid lines 407 and 407 H15 , where we had enough replicates, and tested only for the main effects of parasitoid line, host plant and block.
Likewise, we tested for the main effects of aphid line, parasitoid line and block on a subset of the data including mummies collected from Vicia plants only. We calculated pairwise differences (Tukey HSD) between categories of significant predictors using the package multcomp v1.4.10 (Hothorn et al., 2008).
A model with quasibinomial errors as described above was also applied to analyse the proportion of aphids surviving until the end of the experiment among the non-parasitized aphids (initially exposed aphids minus parasitized aphids). The log-transformed aphid body weight measures were analysed using a linear model and an analysis TA B L E 1 Analysis of deviance

| Parasitism rate and parasitoid emergence
There were highly significant main effects of aphid line, parasitoid line and host plant on parasitism rates, as well as a significant block effect, both in the complete data set and in the data set restricted to H. defensa-infected aphid lines (Table 1). The H. defensa-free aphids and aphids carrying H15 were on average more susceptible to parasitoids than the aphids with the other two strains of H. defensa, and the observed rates of parasitism were mostly higher on Vicia than on the other two plants (Figure 1). The main effect of parasitoid line largely reflects the low parasitism success of line 09-369, even on aphids without H. defensa. In both analyses, parasitism rates were also strongly dependent on the specific combinations of host and parasitoid lines; i.e., there was a highly significant aphid line × parasitoid line interaction, which is uniquely determined by the G × G interaction between L. fabarum and H. defensa in the restricted data set (Table 1). There were no significant effects of the other interaction terms in the full data set, while the parasitoid × plant interaction was marginally significant in the restricted data set. Most notably with regard to the study question, the three-way interaction between aphid line, parasitoid line and host plant was non-significant in both analyses (Table 1). Hence, we observed no evidence for a G × G × E interaction on parasitism rates, indicating that the genetic specificity of the interaction between H. defensa and L. fabarum is not significantly altered by the different host plant environments.
These conclusions remain unchanged when parasitism rates are calculated as the proportion of parasitized aphids among parasitized and surviving adult aphids at the end of the experiment, i.e. excluding aphids that died of reasons other than parasitism (Table S2).
Parasitoid emergence from parasitized aphids differed among treatments. Considering data from aphid lines 407 and 407 H15 separately, parasitoid emergence was in both cases significantly influenced by the host plant, and to a lower extent also by the parasitoid line in aphid line 407 H15 (Figure 2a, b,  (Figure 2c, Table S3).

| Aphid survival and body weight
From all initially exposed aphid nymphs, one part got mummified, one part survived, and one part died before the end of the experiment for reasons other than visible mummification ( Figure S1). Among the non-parasitized aphids, the proportion of surviving aphids varied significantly, explained by a significant main effect of host plant and experimental block (  (Figure 3, Table S4). On average, aphid weight was highest on Beta (0.704 mg ± 0.233 mg SD), followed by Vicia (0.555 mg ±0.161 mg SD) and lowest on Chenopodium (0.400 mg ± 0.110 mg SD).

| DISCUSS ION
Genotype-by-genotype interactions between the parasitoid L. fabarum and the aphid-protective endosymbiont H. defensa have been observed in multiple laboratory experiments (Cayetano and Vorburger, 2015;Schmid et al., 2012) and are therefore assumed to be an important driver of the coevolutionary dynamics in this host-parasitoid system (Hafer & Vorburger, 2019;Hafer-Hahmann & Vorburger, 2020;Kwiatkowski et al., 2012). Here, we investigated such G × G interactions on different host plants, a key environmental variable for aphids and their parasitoids in natural populations. We observed the expected G × G interactions on all host plants, but also a strong main effect of the host plants on overall parasitism rates.
In contrast, we saw no significant G × G × E interaction, suggesting that environmental heterogeneity generated by the availability of different host plants in the field does not reduce the hierarchy of G × G interactions between L. fabarum and H. defensa. Taken together, these findings imply that host plant variation across space and time can indeed create a mosaic of varying selection strength (Thompson, 2005), but without changing the specificity of reciprocal selection between symbiont-protected hosts and parasitoids.
Coevolutionary dynamics may proceed at different pace in different host plant environments, but the environmental variation is unlikely to change the direction of selection (Wolinska & King, 2009 Figure 1). This may not be surprising considering the mechanistic basis underlying H. defensa-conferred resistance, which is related to the presence of a bacteriophage, called APSE, within the bacterial genome (Oliver et al., 2009). Different APSE types carry specific toxin cassettes, which encode for different putative toxins, likely responsible for variation in the protective phenotype among H. defensa strains (Degnan & Moran, 2008a;Moran, Degnan, et al., 2005;Oliver et al., 2009;Oliver & Higashi, 2019). The three strains used here also represent clearly distinct genotypes (Kaech et al., 2021). The phage toxins are assumed to kill susceptible parasitoids at an early stage, that is as eggs or early larvae, but they may also have later-acting effects when parasitoids manage to complete development despite the presence of H. defensa, such as reduced adult weight or delayed emergence of parasitoids (Dennis et al., 2017;Schmid et al., 2012). Such a lateacting detrimental effect may explain the low parasitoid emergence rate from mummies of the aphid line 407 H76 (Figure 2, Table S2), which was already observed in an earlier experiment .
Infection with a toxin-producing symbiont can also be associated with costs to the host, and indeed, the frequency of H. defensacarrying aphids tends to decline within aphid populations that are not under selection by parasitoids (Dykstra et al., 2014;Hafer-Hahmann & Vorburger, 2020;Oliver et al., 2008). Here, we measured adult weight as a rough proxy for aphid performance and found that  Figure S1). Moreover, differences in the reproductive success of parasitoids due to lower parasitism rates on Beta and Chenopodium were amplified further by variation in emergence rates, which were also lowest on Beta and intermediate on Chenopodium Table S2). In summary, Vicia was the most favourable host plant for aphids as well as parasitoids, while Beta represented a comparatively adverse environment for both antagonists. We could thus expect that the strength of reciprocal selection between hosts and parasitoids is higher on relatively benign host plants such as Vicia.
How do differences in host plant quality for both antagonists come about mechanistically? For the parasitoids, differences in parasitism success may be related to variation in plant structure affecting how efficiently aphids are attacked (Grevstad & Klepetka, 1992;Kareiva & Sahakian, 1990), or to variation in host quality (Pan et al., 2020). The entire parasitoid development takes places within the aphid's body; hence, the more vital and well-fed the aphids are, and the more resources may be available for the parasitoid. A resource deficit compared with Vicia-feeding aphids is a possible explanation for the low parasitism success on low-weight aphids from Chenopodium, but another explanation must apply to the frequent parasitoid failure on aphids from Beta, which were even heavier on average than the aphids from Vicia. Beta leaves may contain high concentrations of oxalates (Baker & Eden, 1954), which can have negative effects on aphids (Massonié, 1980) and could thus be responsible for the lower survival in our experiment, but of course many other reasons are conceivable. If the lower vitality of Betafeeding aphids came from the uptake of some toxic plant compound, this may also have hampered the parasitoids' development (Turlings & Benrey, 1998). The low aphid survival on Beta in our experiment was surprising, though, since A. fabae is known as a severe pest in sugar beet cultures (Blackman & Eastop, 2000), and we indeed regularly observe heavy infestations of sugar beet during fieldwork.
F I G U R E 3 Adult weight of the mothers of the experimental aphid generation, in milligrams. Aphids on Beta had the highest weight on average, aphids on Chenopodium the lowest. Boxplot hinges correspond to the 1st and 3rd quartiles, the whiskers extend to a length of 1.5 times the inter-quartile range.
x-axis: endosymbiotic H. defensa strain associated with the aphid clone 407 However, the species Beta vulgaris unites a multitude of crop varieties which may differ in their quality as host plants for A. fabae. The green chard variety we used in this experiment may be a less favourable host plant than the commonly grown varieties of sugar beet.
Whatever the precise reasons, Vicia, Chenopodium and Beta represented very different environments for the aphids and the parasitoids, and still, the genetic interactions tested within remained virtually unaffected. This is not self-evident, as several studies examining host-parasite interactions under different environmental conditions have reported significant G × G × E effects (Bryner & Rigling, 2011;Piculell et al., 2008;Sadd, 2011;Tétard-Jones et al., 2007;Wendling et al., 2017;Zouache et al., 2014). Explicit reports of the lack of a G × G × E interaction in host-parasite systems are scarce to our knowledge (but see Cisarovsky et al. 2012 andVorburger 2013b). To some extent, this may reflect a publication bias against reporting negative results (Csada et al., 1996). We argue that the absence of a significant G × G × E interaction is relevant here because it corroborates the importance we may attribute to the G × G effects. The more robust G × G interactions are to environmental variability, the more pervasive they will be in natural populations, and thus, the more likely they explain fundamental evolutionary phenomena, such as the maintenance of genotypic diversity (Hafer & Vorburger, 2019;Judson, 1995). G × G interactions between H. defensa and L. fabarum have been shown to persist over a range of average temperatures and, as we newly report here, aphid host plants. This consolidates the role of the defensive symbiont H. defensa as a key mediator of coevolution between aphids and parasitoids, also in a heterogeneous environment.

ACK N OWLED G EM ENT
We would like to thank Paula Rodriguez for the diligent long-term maintenance of the insect lines used in this study. This research was funded by the Swiss National Science Foundation (grant no.

31003A_181969 to CV). Open Access Funding provided by Lib4RI
Library for the Research Institutes within the ETH Domain Eawag Empa PSI and WSL.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

AUTH O R CO NTR I B UTI O N S
EG and CV designed the study. EG conducted the experiment. Both authors contributed to data analysis and interpretation. EG wrote the first draft of the manuscript, which was edited and revised by both authors.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.13953.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data set generated in this study is available at Dryad Digital