Positive association between the diversity of symbionts and parasitoids of aphids in field populations

Parasites and pathogens are crucial in shaping immune systems. Many animals and especially insects have outsourced part of their immune function to protective symbionts. There is good evidence that, akin to immune systems, parasites shape the occurrence and diversity of protective symbionts and that likewise, protective symbionts can shape the occurrence and diversity of parasites. Such a relationship should result in a correlation between symbiont and parasite diversity in nature. Aphids are well known for possessing symbionts that provide specific and effective protection against parasitoid wasps. We compared symbiont and parasitoid diversity across multiple populations of different aphid species of the genus Aphis and their parasitoid wasps. The diversity of protective symbionts and parasitoids was indeed positively associated. Even though this association was very noisy, it is in line with the hypothesis that parasitoids and symbionts promote each other’s diversity.


INTRODUCTION
Many animals, especially insects, harbor protective symbionts. Such symbionts can provide effective and specific protection against parasites and pathogens , Fl orez et al. 2015, Vorburger and Perlman 2018; they take over part of the immune system's role. The immune system is one of the most diverse components in any living organism and this diversity is a response to the multitude of different parasites and pathogens it has to cope with (e.g., Apanius et al. 1997, Ghosh et al. 2011. Protective symbionts should be subject to similar selection pressures as a host's own immune system, and therefore, it seems likely that parasites and pathogens may promote patterns of diversity in protective symbionts and vice versa (Hafer and Vorburger 2019). Aphids are particularly well known for their protective symbionts that provide resistance against parasitoids and fungal pathogens (reviewed by Oliver et al. 2014, Vorburger 2014, Guo et al. 2017, Vorburger and Perlman 2018. To date, at least nine species of heritable facultative symbionts have been found in aphids, seven of which have been shown to confer some protection against entomopathogenic fungi or parasitoids (reviewed by Guo et al. 2017). They tend to occur at intermediate frequencies in natural populations, vary within and between aphids species, and can co-occur within the same population or even individual (Russell et al. 2013, Vorburger andRouchet 2016).
The presence or absence of symbionts can drive the abundance of parasitoids and the composition of parasitoid communities both in the laboratory (Sanders et al. 2016, Kraft et al. 2017, K€ ach et al. 2018) and in the field , Rothacher et al. 2016. Protective symbionts impose selection on parasitoids, forcing them to evolve to overcome resistance (Dion et al. 2011, Rouchet and Vorburger 2014, Dennis et al. 2017. By exposing aphids that only differed in whether or not they harbored the protective symbiont Hamiltonella defensa in the field, Rothacher et al. (2016) could show that H. defensa altered the incidence and species composition of naturally colonizing parasitoids. In turn, parasitoids are well known to drive the frequency of protective symbionts (e.g., K€ ach et al. 2018, Ives et al. 2020, Rossbacher and Vorburger 2020. In a recent laboratory experiment exposing aphid populations containing multiple strains of H. defensa to different combinations of asexual lines of the parasitoid wasp Lysiphlebus fabarum, we could show that parasitoid diversity can be crucial in maintaining symbiont diversity (Hafer-Hahmann and Vorburger 2020).
In combination, these findings suggest that parasitoids and protective symbionts of aphids may reciprocally promote diversity (see Hafer and Vorburger 2019), which should result in corresponding patterns observable in nature. To obtain data on parasitoid and symbiont species richness and diversity, we collected aphids from the genus Aphis and their parasitoid wasps from different populations in Switzerland and neighboring countries and combined these samples with samples collected previously for other purposes (Vorburger andRouchet 2016, Vorburger et al. 2017). In line with our expectations, we found that symbiont and parasitoid diversity were indeed positively correlated.

Study system
In this study, we focus on four different aphid species and two subspecies of the genus Aphis: Aphis fabae fabae, Aphis fabae cirsiiacanthoides, Aphis hederae, Aphis ruborum, and Aphis urticata. Each of these species and subspecies specializes on a different host plant or multiple closely related plants, albeit both subspecies of A. fabae share an overwintering host. Aphids in general, and A. fabae in particular, represent a well-known model system for studying host-parasite-symbiont interactions (reviewed by, e.g., Oliver et al. 2014, Vorburger 2014, both under laboratory and under field conditions. Aphis spp. harbors one obligate and several species of facultative bacterial symbionts, some of which offer protection against parasitoids (reviewed by Guo et al. 2017). Facultative endosymbionts are usually transmitted vertically, albeit horizontal transmission may occur occasionally (e.g., Vorburger 2012, Rock et al. 2018). The same symbiont species are usually shared among different species of aphids, and a single aphid can harbor multiple symbiont species , Guo et al. 2017). Frequencies of different symbionts can differ strongly between different species and subspecies of Aphis (Vorburger andRouchet 2016, Vorburger et al. 2017) and vary between populations of the same species (Vorburger and Rouchet 2016). Parasitoids, too, are often shared among different species of aphids albeit the degree of specialization differs strongly between parasitoid species (Raymond et al. 2016).

Field sampling
We collected aphids and parasitoids in June and July 2018 and 2019 in different locations in Switzerland and one location in Germany (see Appendix S1: Figure S1 for map). Five aphid species and subspecies were sampled from six different plants: Aphis fabae fabae from goosefoot (Chenopodium album); Aphis fabae cirsiiacanthoides from thistles (Cirsium arvense and Cirsium vulgare); Aphis hederae from ivy (Hedera helix); Aphis ruborum from black berry (Rubus fruticosus); and Aphis urticata from nettle (Urtica dioica). We followed the same sampling strategy as in Vorburger and Rouchet (2016) and Vorburger et al. (2017), which allowed us to subsequently pool our data with data on aphid symbionts and parasitoids from these studies, which were based on a large collection from 2009. Shortly, to sample aphids, we collected a single, usually adult, aphid from each colony that we found. Usually, an aphid colony is composed of the offspring of a single foundress and hence each aphid within a colony should harbor the same symbiont repertoire. We made sure that we only collected aphids that were at least five meters apart in order to avoid collecting aphids from the same clone, which also resulted in collecting aphids from different plant individuals. Collected v www.esajournals.org aphids were stored at À20°C prior to DNA extraction for symbiont identification (see Parasitoid and symbiont identification). To obtain parasitoids, we collected small plant parts infested with aphid colonies that were visibly parasitized, as indicated by the presence of mummies (dead aphids containing the cocoon of a pupating parasitoid). We brought these colonies to the laboratory and allowed parasitoids to hatch for subsequent identification. For a complete list of samples obtained for each symbiont and parasitoid species, please refer to Appendix S1: Table S1. Our data did not allow us to quantify aphid densities, parasitoid attack rates, or abundance, since the sampling was opportunistic and the focus was on diversity only.

Parasitoid and symbiont identification
Parasitoids were identified using several keys (Japoshvili and Abrantes 2006, Kavallieratos et al. 2013, Hull e et al. 2020. We identified all primary parasitoids to the species level. For a few individuals of Binodoxys and all of Aphelinus collected in 2009, information on species was unavailable. For Binodoxys, we assumed these individuals to belong to the more frequent species in this location and on this host, and for Aphelinus, we assumed all individuals to belong to A. chaonia since this was the only species we encountered in 2018 and 2019. All secondary parasitoids were identified to the genus level. We assume that if parasitoid and symbiont diversity influence each other, such an effect should predominantly occur in primary parasitoids which are directly affected by possible symbionts. However, it is possible that some effects are carried over to secondary parasitoids or that they are influenced by primary parasitoid composition . To identify symbionts, we extracted aphid DNA using high salt extractions following the protocol in Sunnucks and Hales (1994), which was subsequently used for diagnostic PCRs (see Ferrari et al. 2012). We tested for five facultative symbionts: Hamiltonella defensa, Regiella insecticola, Serratia symbiotica, Fukatsuia symbiotica (formerly referred to as X-type/PAXS), and Rickettsia. Four of these, H. defensa (Oliver et al. 2003, Schmid et al. 2012, Asplen et al. 2014), R. insecticola (Vorburger et al. 2010), S. symbiotica (Oliver et al. 2003), and F. symbiotica (Heyworth andFerrari 2015, Donald et al. 2016) have been shown to provide some level of protection against at least some parasitoid species, although the protection may vary among strains of the same species. Additionally, we tested for Buchnera aphidicola, which is an obligate symbiont of aphids and hence served as a control for the success of DNA extraction. Samples that were negative for B. aphidicola were discarded. Please refer to Appendix S1: Table S2 for a full list of primers used. We did not test for three other facultative endosymbionts of aphids, Rickettsiella, Spiroplasma, and Arsenophonus, because we lack data on those from the 2009 collection. While Rickettsiella and Spiroplasma are either absent or exceedingly rare in all study species, Arsenophonus does occur at a substantial frequency (~60%) in one species, A. ruborum (Brechb€ uhler 2019). Hence, we lack one important component of the symbiont community of this species.

Statistical analysis
We used R, version 3.6.1 (R Core Team 2019), to conduct statistical analysis. For all analyses, we included only populations (i.e., samples of the same aphid species, site, and year) for which we had at least six aphid samples, irrespective of the number of parasitoid samples obtained (113 out of 127 populations sampled; average sample sizes [mean AE SD]: 18.43 AE 5.83 aphids and 6.92 AE 5.76 parasitoid samples with on average a total of 84.96 AE 93.87 primary parasitoid individuals per population). This choice represented a compromise between excluding populations for which sampling effort was extremely low and keeping as many populations as possible in the data set. For symbionts, we mostly focused on those symbionts for which at least some protective function against parasitoids has previously been reported (H. defensa, R. insecticola, S. symbiotica, and F. symbiotica). We looked at diversity using two different measures, species richness and species evenness (Shannon diversity) which we calculated using the R package vegan (Oksanen et al. 2019). Since some of our sample sizes were not high enough to reach saturation and hence sample size may have explained some of the variation in these estimates based on raw numbers, we additionally calculated rarefied diversity estimates. To do so, we only used v www.esajournals.org populations in which we identified a minimum of six symbionts and for which we had at least 6 parasitoid samples with 6 or more different (i.e., different species or in different samples) primary parasitoid occurrences and more than ten primary parasitoid individuals in total (42 populations; average sample sizes [mean AE SD]: 19.14 AE 5.32 aphids and 10.95 AE 5.56 parasitoid samples with on average a total of 134.29 AE 106.12 primary parasitoid individuals per population). Again this represents a trade-off between ensuring reasonable sampling effort for each population and keeping as many populations as possible in the data set. To obtain rarefied diversity estimates for Hill numbers 0 (i.e., species richness) and 1 (i.e., Shannon diversity), we used estimateD in iNEXT (Chao et al. 2014) and specified a sample size of 6 due to our cutoff of 6 samples. For Shannon index and rarefaction, we focused on the number of samples (i.e., aphid colonies) in which each parasitoid occurred since the number of aphids, and hence, mummies can differ greatly between colonies and would have a large impact on the number of individuals produced in each colony. For each response variable, we built a linear mixed model (Bates et al. 2015) using the appropriate measure of symbiont diversity as response and the corresponding parasitoid diversity as predictor. To analyze the unrarefied estimates of species richness and diversity, we additionally included the number of aphid samples as a covariate. We chose to use linear mixed models rather than simpler correlation tests to be able to account for aphid host species, site, and year as random effects. Unlike a simple correlation, these models imply a direction, that is, we tested for an effect of parasitoid diversity on symbiont diversity. We chose this direction since we could show previously that under laboratory conditions, parasitoid diversity can indeed affect symbiont diversity (Hafer-Hahmann and Vorburger 2020). However, we are aware that due to the correlative nature of our findings, we cannot rule out that any significant effect of parasitoid diversity we observe may rather be an effect of symbiont diversity on parasitoid diversity or be caused by different factors jointly influencing both symbiont and parasitoid diversity. We used TransformTukey from package rcompanion (Mangiafico 2019) on each response variable prior to analysis to comply with model assumptions. In order to test whether environmental factors affected symbiont diversity (annual mean temperatures and annual precipitation), we included these environmental variables as fixed factors in separate linear mixed effect models. To analyze a correlation between parasitoid diversity and within aphid symbiont diversity, we calculated the mean number of different symbionts within each aphid and used linear mixed models as described above with each of our four estimates of parasitoid diversity as fixed factor. Each model was then followed with a type III analysis of variance (ANOVA) using lmerTest to obtain p values with Satterthwaite's degrees of freedom method (Kuznetsova et al. 2017). Plots were created using ggplot2 (Wickham 2016).

RESULTS
Primary parasitoid diversity had a significant effect on protective symbiont (H. defensa, R. insecticola, S. symbiotica, and F. symbiotica) diversity in our models (species richness F 1,108 = 9.00, P = 0.0034, with a significant effect of the number of aphids sampled F 1,108 = 5.35, P = 0.0226; Shannon index F 1,109 = 15.24, P = 0.0002, without any significant effect of the number of aphids sampled F 1,103 = 0.36, P = 0.5474), albeit this relationship was very noisy ( Figure 1A,B). These findings remained robust when we used rarefied diversity estimates instead (species number F 1,39 = 8.70, P = 0.0054; Shannon index F 1,39 = 8.59, P = 0.0056; Figure 1C,D). This positive association with protective symbiont diversity also held true when we considered all parasitoids, that is, included secondary parasitoid diversity (Appendix S1: Table S4A). When we used symbiont diversity calculated for all symbionts investigated in this study (H. defensa, R. insecticola, S. symbiotica, F. symbiotica, and Rickettsia), the association between primary parasitoid and symbiont diversity remained similar, albeit less pronounced, when using raw diversity estimates, but disappeared when we used rarefied diversity estimates instead (Appendix S1: Figure S2, Table S4B). The latter was arguably due to an increase of rarefied symbiont diversity in one of the aphid species, A. hederae, which was host to a comparatively low number of different parasitoid species (Appendix S1: Figure S2 C,D).  We found no significant associations between precipitation or annual temperatures and symbiont diversity (Appendix S1: Table S5). The diversity of symbionts within aphids (i.e., species richness of different symbionts within each aphid individual), if anything, tended to also be positively associated with primary parasitoid diversity (Appendix S1: Figure S3, Table S6).

DISCUSSION
The diversity of protective symbionts and parasitoid wasps was correlated across populations of different aphid species of the genus Aphis. This held true not only for primary parasitoids, which we know to interact with symbionts, but also when we included secondary parasitoids. Even though the association was very noisy, it is in line with the hypothesis that protective symbionts and parasites may drive each other's diversity (reviewed by Hafer and Vorburger 2019). This hypothesis is supported by experimental evidence showing that a more diverse population of parasitoids will maintain a higher diversity of symbiont strains in aphids (Hafer-Hahmann andVorburger 2020, Rossbacher and. Our findings complement these experiments and indicate that the processes found to shape the interaction of symbionts and parasites under more artificial conditions may indeed play a role in shaping natural communities. However, the correlative nature of our study does not allow us to infer the directionality of the effect-protective symbionts may also have bottom-up effects on parasitoid diversity (Rothacher et al. 2016)nor does it allow to rule out alternative explanations for our findings. Biotic or abiotic environmental factors jointly influencing both parasite and symbiont diversity could also generate the observed association, but this would not be mutually exclusive with the explanation proposed above.
Interestingly, the association we observed was much clearer if we excluded Rickettsia which, unlike the other symbionts we used, has not previously been observed to protect against parasitoids. This seems to support our hypothesis which should apply exclusively to symbionts protective against parasitoids. However, the evidence for the protective ability of the symbionts we analyzed also showed important variation among strains (e.g., Hansen et al. 2012). We can therefore not exclude that a protective function may be found for certain strains of Rickettsia in the future, which could invalidate our distinction of protective versus non-protective symbionts.
All aphid species we studied are cyclical parthenogens, that is, they reproduce clonally during summer, and their symbionts are predominantly vertically transmitted. Hence, we might expect some association between clonal diversity and the diversity of protective symbionts. Any ecological process promoting clonal diversity in aphid populations could affect symbiont diversity indirectly (although heritable symbionts are of course part of the clonal genotypes). However, the protective effect provided by symbionts is usually much stronger than variation in the aphids' innate resistance to parasitoids (e.g., Oliver et al. 2005, Martinez et al. 2014, Vorburger 2014. Therefore, we expect that symbionts play a larger role than aphid genotypes when it comes to interactions with parasitoids. As mentioned above, protection can also vary greatly between strains of the same symbiont species (e.g., Hansen et al. 2012, Cayetano and. We did not measure strain diversity of symbionts and hence may have left out a functionally important aspect of symbiont diversity. This may help explain some of the noise observed in our data. Our understanding of the naturally occurring strains differs between different symbiont species, but we know that multiple strains can co-occur within a single population (Henry et al. 2013, Brechb€ uhler 2019, Oliver and Higashi 2019. It is plausible that strain-level diversity varies among different symbiont species, and hence, by ignoring strain diversity we may have underestimated (functional) diversity for some symbionts more than for others. Similarly, we also did not look at within parasitoid diversity; at least for some parasitoid species, different genotypes vary in their ability to overcome resistance conferred by specific symbiont strains (Rouchet and Vorburger 2012, Schmid et al. 2012, Cayetano and Vorburger 2013, 2015, Vorburger and Rouchet 2016. Given that we did not consider functional diversity that was likely present within symbiont as well as parasitoid species, it is interesting that a positive v www.esajournals.org association was observed at the coarser, specieslevel diversity we assessed. This should warrant further study both on the finer, within-species diversity level and under more controlled laboratory conditions. Protection against parasitoids is not the only benefit symbionts confer to their hosts. Regiella insecticola (Scarborough et al. 2005), F. symbiotica (Heyworth and Ferrari 2015), and Rickettsia (Łukasik et al. 2013) protect also against entomopathogenic fungi and can do so in a specific manner . Abiotic factors, too, have the potential to affect symbiont composition. Hamiltonella defensa, R. insecticola, and F. symbiotica, for example, have been found to protect against heat shock (e.g., Moran 2006, Heyworth andFerrari 2015), while the costs of S. symbiotica increase following a heat shock (Russell and Moran 2006). Surprisingly, we did not find any significant associations between symbiont diversity and temperature or precipitation. Such associations were detected in previous studies on aphids (Tsuchida et al. 2002, Sep ulveda et al. 2017) and other arthropods (Toju andFukatsu 2011, Zhu et al. 2018). It is possible that our measure of mean annual temperature did not accurately represent the risk of extreme temperatures or that environmental factors did not span a sufficient range across our study sites to identify any clear associations.
Earlier studies investigating protective symbionts in a geographic context have focused on correlations between symbiont prevalence and parasitoid incidence, with results varying between studies, populations, and host plants (Hansen et al. 2007, Smith et al. 2015, Zepeda-Paulo et al. 2016. Recently, Ye et al. (2018) observed associations between two protective symbionts, H. defensa, and R. insecticola in the grain aphid and parasitoid communities. The data we obtained were not suitable to judge the abundance of parasitoids. Rather, we focused on their diversity and did so over a relatively large geographic range. To our knowledge, this presents the first study explicitly looking at the correlation between symbiont and parasite diversity in the field. However, correlations between immune diversity and parasite diversity have previously been observed in other systems, such as in the major histocompatibility complex (MHC) of vertebrates (e.g., Go€ uy De Bellocq et al. 2008, Eizaguirre et al. 2011, Minias et al. 2020, supporting the idea that protective symbionts could be under similar selection as host immune genes (Hafer and Vorburger 2019).
Similar correlations with regard to species diversity are common in ecology, whereby bottom-up and/or top-down effects of diversity at one trophic level are expected to facilitate diversity on another trophic level. Examples for this come from systems as diverse as mammals and parasites (e.g., Nunn et al. 2004, Morand 2015, plants, herbivores, their parasitoids, and mutualists (Mitter et al. 1991, Mailafiya et al. 2010, Cao et al. 2018, including parasitoid host associations (Whitfield 1998). Here we found that similar effects may be at play when it comes to protective symbionts and the parasitoids against which they provide protection.
Assuming that parasitoid diversity indeed played a role in shaping the symbiont diversity in the aphids we studied, one can ask about ecological factors influencing parasitoid diversity. An interesting study by  found that higher plant diversity can enhance population level symbiont diversity in aphids. While this may reflect a direct effect of plant diversity, an indirect effect via parasitoids would be plausible as well. We used aphid species highly specialized on certain plant species and focused mainly on a single host plant for each aphid species, hence ignoring potential differences in symbiont communities between host plants. Such differences are indeed known from other aphid species, such as the pea aphid Acyrthosiphon pisum (Ferrari et al. 2004(Ferrari et al. , 2012, or the cowpea aphid Aphis craccivora (Brady and White 2013). The generally high host plant specificity of aphids entails that more diverse plant communities support more aphid species (e.g., Joshi et al. 2004), and species-level aphid diversity correlates positively with parasitoid diversity (Petermann et al. 2010). We did not characterize the available host plants and their aphid communities in the vicinity of the focal aphid populations we studied, but it is plausible that the surrounding communities influenced the patterns we observed. This would suggest an interesting scenario where diversity begets diversity across multiple levels, namely from plants to aphids, their symbionts, and parasitoids.