A Link Between Communities of Protective Endosymbionts and Parasitoids of the Pea Aphid Revealed in Unmanipulated Agricultural Systems

In the last decade, the influence of microbial symbionts on ecological and physiological traits of their hosts has been increasingly recognized. However, most of these effects have been revealed under laboratory conditions, which oversimplifies the complexity of the factors involved in the dynamics of symbiotic associations in nature. The pea aphid, Acyrthosiphon pisum, forms a complex of plant-adapted biotypes, which strongly differ in the prevalence of their facultative endosymbionts. Some of the facultative endosymbionts of A. pisum have been shown to confer protection against natural enemies, among which Hamiltonella defensa is known to protect its host from parasitoid wasps. Here, we tested under natural conditions whether the endosymbiont communities of different A. pisum biotypes had a protective effect on their hosts and whether endosymbiotic associations and parasitoid communities associated with the pea aphid complex were linked. A space-time monitoring of symbiotic associations, parasitoid pressure and parasitoid communities was carried out in three A. pisum biotypes respectively specialized on Medicago sativa (alfalfa), Pisum sativum (pea), and Trifolium sp. (clover) throughout the whole cropping season. While symbiotic associations, and to a lesser extent, parasitoid communities were stable over time and structured mainly by the A. pisum biotypes, the parasitoid pressure strongly varied during the season and differed among the three biotypes. This suggests a limited influence of parasitoid pressure on the dynamics of facultative endosymbionts at a seasonal scale. However, we found a positive correlation between the α and β diversities of the endosymbiont and parasitoid communities, indicating interactions between these two guilds. Also, we revealed a negative correlation between the prevalence of H. defensa and Fukatsuia symbiotica in co-infection and the intensity of parasitoid pressure in the alfalfa biotype, confirming in field conditions the protective effect of this symbiotic combination.

In the last decade, the influence of microbial symbionts on ecological and physiological traits of their hosts has been increasingly recognized. However, most of these effects have been revealed under laboratory conditions, which oversimplifies the complexity of the factors involved in the dynamics of symbiotic associations in nature. The pea aphid, Acyrthosiphon pisum, forms a complex of plant-adapted biotypes, which strongly differ in the prevalence of their facultative endosymbionts. Some of the facultative endosymbionts of A. pisum have been shown to confer protection against natural enemies, among which Hamiltonella defensa is known to protect its host from parasitoid wasps. Here, we tested under natural conditions whether the endosymbiont communities of different A. pisum biotypes had a protective effect on their hosts and whether endosymbiotic associations and parasitoid communities associated with the pea aphid complex were linked. A space-time monitoring of symbiotic associations, parasitoid pressure and parasitoid communities was carried out in three A. pisum biotypes respectively specialized on Medicago sativa (alfalfa), Pisum sativum (pea), and Trifolium sp. (clover) throughout the whole cropping season. While symbiotic associations, and to a lesser extent, parasitoid communities were stable over time and structured mainly by the A. pisum biotypes, the parasitoid pressure strongly varied during the season and differed among the three biotypes. This suggests a limited influence of parasitoid pressure on the dynamics of facultative endosymbionts at a seasonal scale. However, we found a positive correlation between the α and β diversities of the endosymbiont and parasitoid communities, indicating interactions between these two guilds. Also, we revealed a negative correlation between the prevalence of H. defensa and Fukatsuia symbiotica in co-infection and the intensity of parasitoid pressure in the alfalfa biotype, confirming in field conditions the protective effect of this symbiotic combination.

INTRODUCTION
Antagonistic interactions such as host-parasite or prey-predator are major drivers of ecological and evolutionary processes and affect more globally biodiversity patterns and ecosystem functioning. Antagonistic interactions usually involve arms race between protagonists whereby the enemy evolves new weapons to counteract victim's defensive strategies (Abrams, 2000). Selection of these adaptive responses acts on variation encoded by enemy or victim's genomes and relying on different mechanisms such as toxins released by pathogens into the infected host or immunity host defenses triggered upon parasite attacks. More recently, evidence accumulated that some microbial symbionts hosted by eukaryotes can also actively participate in antagonistic interactions by extending the arsenal of enemies or the defensive strategies of victims (Flórez et al., 2015).
Protective microbial symbionts have been particularly wellstudied in insects with the best examples being the Spiroplasma bacterial endosymbiont which protects a mushroom-feeding fly Drosophila neotestacea from parasitic nematodes (Jaenike et al., 2010) or the Hamiltonella defensa bacterial endosymbiont which confers protection to aphids against hymenopteran parasitoids (Oliver et al., 2003). In both cases, the mechanisms of the symbiont-mediated protection involve the production of microbial toxins active against the parasite (Oliver and Perlman, 2020). The infection with these protective symbionts confers an immediate fitness advantage to the host when enemy pressure exists. However, hosting protective symbionts may come with some fitness costs (Oliver et al., 2008;Simon et al., 2011;Vorburger and Perlman, 2018). This cost/benefit balance has been repeatedly invoked to explain intermediate frequencies of protective symbionts in natural populations of their hosts (Oliver et al., 2014). In addition, these bacterial symbionts being transmitted occasionally through horizontal transfer events, they may rapidly spread across both host populations and species, leading to rapid adaptations and invasion processes (Jaenike et al., 2010;Himler et al., 2011). Finally, these protective symbionts may modify the diversity, structure and function of ecological networks that link their hosts to lower or higher trophic levels, through bottom-up and top-down effects (Mclean and Godfray, 2015;Rothacher et al., 2016;Ye et al., 2018;Mclean, 2019). Most studies on protective symbionts in insects have been carried out in laboratory conditions or in simplified field experiments. Although these studies have been crucial in investigating the costs and benefits in hosting protective symbionts in presence or absence of parasitism, only few assessed the response of host-symbiont associations to varying pressure intensities of natural enemies and by extent, the link between symbiont communities and parasitoid communities in real and complex natural environments (Smith et al., 2015;Ye et al., 2018).
The pea aphid Acyrthosiphon pisum is a good model to study how symbiotic associations are influenced by the parasitism exerted on their hosts and in return, how symbionts affect diversity and community structure of their host's natural enemies. This sap-feeding insect feeds on legumes and forms at least 15 biotypes differentiated by host plant utilization (Peccoud et al., 2009(Peccoud et al., , 2015Ferrari et al., 2012). A large variation in symbiotic associations is found in natural populations of pea aphids with seven heritable facultative endosymbionts being hosted alone or in co-infections in addition to the primary symbiont Buchnera aphidicola (Gauthier et al., 2015;Guyomar et al., 2018). These endosymbionts can have different phenotypic effects on their hosts, including protection from natural enemies. For example, H. defensa, Fukatsuia symbiotica (also named PAXS for Pea Aphid X-type Symbiont), and Serratia symbiotica have been reported as protective against parasitoids (Oliver et al., 2003;Ferrari et al., 2004;Guay et al., 2009;Leclair et al., 2016), while Regiella insecticola, Rickettsiella viridis, Rickettsia sp., Spiroplasma sp., and F. symbiotica confer protection against fungal pathogens (Scarborough et al., 2005;Łukasik et al., 2013;Heyworth and Ferrari, 2015) and R. viridis against predators (Polin et al., 2015). Beyond the phenotype induced by one endosymbiont, other phenotypes are observed when multiple symbionts coexist in a single host (Leclair et al., 2017;Mclean et al., 2018). For example, multiple infections with H. defensa and F. symbiotica can give stronger protection against parasitoids than each endosymbiont alone (Guay et al., 2009;Heyworth and Ferrari, 2015;Leclair et al., 2016), although this may not be always the case (Doremus and Oliver, 2017). This variation in phenotypic effects induced by the different endosymbionts could result from strain differences as showed in earlier works (Oliver et al., 2003;Cayetano et al., 2015;Leclair et al., 2016).
Population cage experiments showed that pea aphids infected with H. defensa increased rapidly in frequency when exposed to Aphidius ervi, the main parasitoid of A. pisum, but declined when parasitoid pressure was removed (Oliver et al., 2008). Recent independent studies exposed laboratory clones carrying or not H. defensa to field conditions in order to analyze the temporal dynamics of protected and unprotected aphid clones when facing a wider range of natural enemies (Hrcek et al., 2016;Rothacher et al., 2016). While parasitism rates were lower in protected clones, this advantage did not translate into higher population growth rate. Overall, these results indicate that although symbiont protection is expressed in field conditions, it is mitigated by multiple factors that may include the rate of parasitism, the range of natural enemies and the costs entailed by harboring the protective symbionts. The work of Smith et al. (2015) was the first to analyze the dynamics of protective endosymbionts in relation with parasitism pressures in natural populations of the pea aphid. Important seasonal shifts in symbiont composition were found, possibly driven by changes in parasitism rates. In particular, the frequency of H. defensa was related with parasitoid pressure while the frequency of R. insecticola was linked with fungal-induced mortality. Additionally, the aforementioned studies highlighted the influence of protective endosymbionts on diversity and structure of parasitoid communities, with potential consequences on food web functioning Mclean, 2019). More recently, in an elegant experimental evolution study (Hafer-Hahmann and Vorburger, 2020), it was discovered that symbiont diversity could be driven by parasitoid diversity, suggesting potential reciprocal influence of these two communities through their shared aphid hosts.
In the present study, we tested under natural conditions whether the endosymbiont communities of different A. pisum biotypes had an actual protective effect on their hosts and whether endosymbiont communities and parasitoid communities associated with the pea aphid complex were linked. To this end, we measured parasitoid pressure on various pea aphid populations, examined whether facultative symbiotic associations responded to parasitoid pressures of varying intensity, and analyzed the relationships between symbiotic and parasitoid communities in terms of both α-diversity (i.e., the species diversity in a given sample) and β-diversity (i.e., the change in diversity of species from one sample to another). In essence, our study is similar to that of Smith et al. (2015) but differs in that we performed field surveys in the native range of A. pisum and its natural enemies [and not in the introduced range as in Smith et al. (2015)]. We therefore expect a wider diversity in aphid-endosymbiont-enemy interactions that may drive different patterns of protective symbiont dynamics. To address these issues, we carried out a field survey, spanning twelve dates across the growing season, and encompassing three biotypes of A. pisum, each specialized to distinct legume crops and harboring or not protective facultative endosymbionts. For each date and crop, we estimated parasitoid pressure exerted on pea aphids using two proxies, determined the structure of the parasitoid community, recorded the frequencies of each facultative endosymbiont in the pea aphid populations and their distribution among host individuals (i.e., symbiont communities), analyzed how the prevalence of symbiont communities influenced parasitism rate within aphid populations, and measured how diversities of both parasitoid and facultative symbiont communities were related. Finally, because parasitoids may also drive symbiont diversity at the population level (Hafer and Vorburger, 2019), we examined in a subsample of our aphid collections whether variation in the strain of the protective endosymbiont H. defensa was related to the intensity of parasitism pressures.

Sampling Sites
Pea aphid individuals were collected in crop fields of alfalfa (Medicago sativa), clover (Trifolium sp.), and pea (Pisum sativum) in Western France (25 km around Rennes, Brittany) between May and October 2014 during the parthenogenetic phase of the aphid life cycle. Locations of the sampled fields and their coding are given in Supplementary Figure 1 and Supplementary  Table 1. When possible, three distant fields per legume species were sampled every two weeks starting on May 7 and ending on October 8, 2014, resulting in 12 sampling dates (Supplementary Table 2). As in Brittany pea is harvested by the end of July, no sample was available for this legume in August. By chance, some pea regrowth colonized by aphids were found at beginning of September and were sampled until October 8. While we aimed to monitor the same crop fields during our survey, some fields were sometimes not available because of either the absence of insects or the crop harvesting (see Supplementary  Table 1 for details). From here, one sampling refers to an insects' collection made at one date in one field of a given legume species. Overall, we made a total of 34, 31, and 19 samplings for alfalfa, clover, and pea, respectively. At each sampling, alive pea aphid adults and nymphs were collected using beat sampling from plants separated by approximately 5 m to cover a larger spatial/genetic diversity at the field scale. Aphid mummies (i.e., a dead aphid containing an immature parasitoid) were also collected by view in order to extract parasitoid specimens from field.

Estimates of Parasitism Pressure
To estimate the parasitism pressure exerted on aphid populations, we used two different proxies. For the first one, we collected on average 80 adult aphids at each field sampling (see Supplementary Table 2 for a precise number of sampled adults per field) that were then installed by batch of ten on plants of broad bean (Vicia faba), which is used as a host plant by the different pea aphid biotypes (Peccoud and Simon, 2010), under laboratory conditions (20 • C and 16hD:8hN regime). Ten days after the sampling, the rate of parasitism was assessed by calculating the ratio between the number of aphid mummies obtained and the total number of installed aphids for each field per date sample. Overall, 3,074 aphid adults from alfalfa, 1,830 from clover and 1,637 from pea were followed for their parasitism status. A second proxy of parasitism pressure was considered in order to reduce the potential estimate bias caused by aphid resistance against parasitoids. For this second proxy, we sampled the field canopy by using a leaf vacuum cleaner during five seconds at three different sites randomly chosen in a given field. The three canopy samples were kept in a plastic bag separately and stored at −20 • C. In the laboratory, pea aphids (all instars combined) were counted and parasitoids were numbered and preserved in 70% ethanol for species identification. The second proxy of parasitism pressure consisted on calculating the ratio between the number of parasitoids belonging to the guild using pea aphids as hosts (Starý, 2006) and the number of pea aphids.

Aphid Parasitoid Communities
To estimate the diversity and structure of parasitoid communities, we considered parasitoid specimens emerging from pea aphid mummies collected from the fields or obtained in laboratory. Those parasitoid specimens were determined at species level (when possible) based on morphological criteria (Müller et al., 1999;Starý, 2006). The present study excluded all hyperparasitoids from the analyses (258 specimens were found). Overall, we determined 744 primary parasitoids from alfalfa, 794 from clover, and 373 from pea, respectively. From all these specimens, six parasitoid taxa were identified: A. ervi, Aphidius eadyi, Aphidius avenae, Praon volucre, Praon barbatum, and Aphelinus sp.

Diversity and Relative Abundance of Facultative Endosymbionts in the Pea Aphid
For each sampling, 60 pea aphid nymphs were collected in order to determine the symbiotic composition in the sampled aphid population. Young aphid nymphs (second-third instars) were chosen instead of older stages to reduce the bias induced by differential mortality caused by parasitoids. After sampling, the aphid nymphs were preserved in 95% ethanol. A random subsample of about 30 aphids per sampling was then inspected for their bacterial symbiont composition (see Supplementary Table 2 for precise number of aphids inspected per field). A total of 993 of aphids from alfalfa, 831 from clover and 509 from pea were screened for the presence of the seven facultative endosymbionts reported in the pea aphid: S. symbiotica, H. defensa, R. insecticola, F. symbiotica (or "Pea Aphid X-type Symbiont"), R. viridis, Rickettsia sp., and Spiroplasma sp. Their presence was detected using speciesspecific PCR primers according to Peccoud et al. (2014) and the detection of the obligatory bacterial endosymbiont B. aphidicola was used as a control for DNA extraction. Note that we also searched for Wolbachia because of some reports for the pea aphid in the literature (Russell et al., 2013;Wang et al., 2014;Gauthier et al., 2015) but we did not detect any positive aphids. As in the past, we have conducted extensive surveys of the microbial communities associated with the pea aphid, using either targeted or whole-genome metagenomics approaches (Gauthier et al., 2015;Guyomar et al., 2018), which have confirmed the presence of seven facultative endosymbionts considered in our paper, we are confident that we have surveyed the main fraction of the pea aphid symbionts. For each sampling, we thus estimated the prevalence of these seven facultative heritable endosymbionts and their distribution among aphid individuals. The endosymbiont infection status of each aphid individual (i.e., identity and number of endosymbionts) was then determined.

Seasonal Variation in Hamiltonella Strains
In order to monitor possible changes in H. defensa strain diversity throughout the cropping season that could result from a variation in the intensity of parasitoid pressure, a subsample of 72 pea aphid nymphs from the alfalfa biotype and infected with H. defensa singly (30 aphids) or in coinfection with F. symbiotica (42 aphids) was selected among the samplings done before (36 aphids) and after the peak (36 aphids) of parasitoids' activity. All aphid individuals were genotyped at seven microsatellite markers following (Peccoud et al., 2008) in order to distinguish aphid clones based on their multilocus genotype. As the parasitism protection due to H. defensa infection is associated with the presence of a bacteriophage (APSE), which encodes toxins potentially responsible for parasitoid development arrest (Degnan and Moran, 2008), phage presence was checked by PCR using the primer pair P3 (forward) 5 -TCGGGCGTAGTGTTAATGAC-3 (reverse) 5 -TTCCATAGCGGAATCAAAGG-3 and P51 (forward) 5 -AG GTGCGATTACCCTGTTTG-3 (reverse) 5 -GATAAAACATCG CCGTTTGC-3 (Mclean and Godfray, 2015). A multilocus sequence-typing (MLST) was then performed for the characterization of H. defensa strains with partial DNA sequences of housekeeping genes accD and murE (Henry et al., 2013). Fragments were amplified by PCR using H. defensaspecific primers and cycling conditions described in Henry et al. (2013). Amplicons were sent to Genoscreen for Sanger sequencing. Sequences obtained were cleaned and aligned using Geneious R v.7.1.5 (Kearse et al., 2012). For each sample and each gene (accD and murE) sequences were used to build a phylogenetic tree using the Neighbor Joining method (Tamura-Nei distance). Bootstrap values were computed for each branch node (N = 1,000).

Parasitism Proxies and Temporal Variation in Parasitism Pressure
The first analysis consisted in assessing the quality of parasitism pressure estimates by calculating the Spearman correlation coefficient between the two proxies. Then, we analyzed the temporal dynamics of parasitoid activity in the three pea aphid biotypes. From the latter analysis, we found three distinct parasitism periods arbitrarily defined as "pre-parasitism, " "maximal parasitism, " and "post-parasitism." To test whether the overall parasitism rate differed between the three biotypes, we used Generalized Linear Mixed Model (GLMM) (glmer function of the lme4 package (Bates et al., 2015) with a binomial error distribution (logit link function). Both the field ID and session number were fit as random factor to include data substructure.

Assessment of α-Diversity of Parasitoid Communities and Aphid Endosymbiont Infection Statuses
For each combination of parasitism period and pea aphid biotype, the occurrence and co-occurrence of parasitoid species and facultative endosymbionts were visualized using the R package Mondrian (Siberchicot et al., 2016). For each combination biotype/session (i.e., we pooled the data obtained in the fields of a given legume crop during one sampling session), the α-diversity of parasitoid communities was estimated using the Shannon Index that accounts for both abundance and evenness of the species. Also, we calculated the αdiversity of the aphid endosymbiont infection statuses using the Shannon Index by considering the relative abundance of each infection status found in a given biotype/session combination. By using a General Linear Models, we tested whether the parasitism period (i.e., three levels factor: "preparasitism, " "maximal parasitism, " and "post-parasitism") or the pea aphid biotype (i.e., three levels factor: "alfalfa, " "clover, " and "pea" biotypes) affected the Shannon index of parasitoid communities or endosymbiont infection statuses of aphids. Finally, to test whether α-diversities of both parasitoid communities and aphid endosymbiont infection statuses were linked, a Pearson correlation coefficient between their Shannon indexes was calculated.

Assessment of β-Diversity of Parasitoid Community and Aphid Endosymbiont Infection Statuses
We quantified the dissimilarity between the parasitoid communities or the endosymbiont infection statuses of aphids (β-diversity) between all pea aphid samples using the Bray-Curtis distance. The Bray-Curtis dissimilarities between all pairwise combinations of pea aphid samples were calculated, ordinated following a non-metric multidimensional scaling (nMDS) and represented on a scatter graph where the position of an endosymbiont status or a parasitoid community depended on its distance from all other points in the analysis. Then, the effects of the pea aphid biotype (i.e., three levels factor: "alfalfa", "clover, " and "pea" biotypes) and the parasitism period (i.e., three levels factor: "preparasitism, " "maximal parasitism, " and "post-parasitism") on endosymbiont infection statuses of aphids or parasitoid community dissimilarity were tested by performing a permutational multivariate analysis of variance (PERMANOVA) on the Bray-Curtis dissimilarity matrix. Pairwise comparisons between levels of factors were performed using pairwise Adonis tests with Bonferroni corrections (Martinez-Arbizu, 2017). Finally, the correlation between the two dissimilarity matrices (β-diversity of endosymbiont infection statuses of aphids vs β-diversity of parasitoid communities) was calculated by using a Mantel test. The Mantel test, nMDS, and PERMANOVA were implemented using the R package vegan (Oksanen et al., 2019), and the graphical representation of the nMDS was generated using the R package ggplot2 (Wickham, 2016).

Prevalent Symbiotic Associations and Parasitism Rates
To search for a link between symbiont associations and parasitoid pressures, the parasitism rate due to a dominant parasitoid species or all parasitoid species was tested against the frequency of the symbiotic associations in each biotype using General Linear Mixed Models (LMM) (lme function of the nlme package; Pinheiro et al., 2020). Note that we did not test all the correlations as some infection states had too few observations. As some fields were considered several times during our survey, the field ID was fit as a random factor to include data dependency. Before each statistical modeling, we checked linearity assumption between the response and the explanatory covariate. In case of linearity departure, we used Generalized Additive Mixed Models (GAMM) (gamm function of the mgcv package; Wood, 2017). Model assumptions were verified by plotting residuals versus fitted values for each model.

Parasitism and Variation in H. defensa Strain Diversity
To analyze the effect of parasitism pressure on the H. defensa strain diversity, we compared the frequencies of each H. defensa strain found before and after the parasitism peak in alfalfa fields using a χ 2 test.

Parasitism Proxies and Temporal Variation in Parasitism Pressure
The two proxies of parasitism pressure were positively correlated (Spearman correlation: rho = 0.452, p < 0.001) and showed similar dynamics during the survey (Figure 1): a low parasitism (<0.3) early in the season, then a maximal parasitism in June FIGURE 1 | Temporal variation in parasitism pressure in alfalfa, clover, and pea crops estimated using two proxies: parasitism rate (black dots, -) and the ratio of the number of parasitoids/pea aphid individuals (white dots, -). The three shaded areas within each graph correspond to the three parasitism periods arbitrarily defined: "pre-parasitism," "maximum parasitism," and "post-parasitism" periods. Error bars represent standard errors. and July (up to 0.5) and a decrease below 0.2 from late July to October. The same temporal dynamics of parasitism activity was found in the three legume crops and from these results, three distinct periods were arbitrarily defined as "preparasitism, " "maximal parasitism" (0.45-0.65 parasitism rate), and "post-parasitism" (Figure 1). Overall, the rates of parasitism were the highest in clover fields (30.3%), intermediate in pea (22.4%), and the lowest in alfalfa (15.9%) (GLMM: χ 2 = 10.26, df = 2, p < 0.005).

α-Diversity in Parasitoid Communities and Endosymbiont Infection Statuses of Aphids
Overall, A. ervi was the dominant parasitoid species in both alfalfa and clover crops (respectively 79.23 and 94.29% of emerging parasitoids over all the season) and the second most abundant species in pea crops (36.44%), after A. eadyi (45.48%) and before A. avenae (15.16%) (Figure 2). Both P. barbatum and Aphelinus sp. were exclusive to the clover and alfalfa biotypes, with P. barbatum accounting for 7.34% of parasitoids in alfalfa crops. The Shannon index in the parasitoid communities varied according to the pea aphid biotype and the parasitism period (LM: interaction term, F = 4.00, df = 4, p = 0.013). Overall, the α-diversity in parasitoid communities was the lowest in the clover fields, intermediate in alfalfa and the highest in pea. Although the dominant parasitoid species remained the same during the season in alfalfa and clover biotypes, species richness of their parasitoid communities increased throughout the season (Figure 3). In pea crops, A. eadyi increased in frequency, reaching 82% of parasitoids in post-parasitism period, replacing A. ervi as dominant species in the guild and leading to a decline in the Shannon index at this period.
Overall, 94.7% of pea aphids (96.5% for aphids from alfalfa, 94.7% from clover and 91.3% from pea) were infected by at least one facultative bacterial endosymbiont (Figure 4). Seven facultative endosymbiont species were detected with prevalence varying strongly between A. pisum biotypes. On average aphid individuals harbored 1.8, 1.1, and 1.4 facultative endosymbionts in alfalfa, clover and pea crops, respectively. In alfalfa, coinfection of H. defensa with F. symbiotica was the most frequent symbiotic association across the field survey (accounting for 43.70% of aphids) followed by H. defensa in monoinfection (12.28%). Two other associations involving H. defensa were also noted in alfalfa fields: coinfection with R. insecticola (8.96%) and triple infection with R. viridis and F. symbiotica (7.25%). Aphids free of secondary endosymbionts in alfalfa represented only 3.5%. In clover fields, R. insecticola in monoinfection represented 74.48% of surveyed aphids while individuals free of any facultative endosymbiont represented only 5.29%. In pea fields, infection with S. symbiotica singly (37.25%) or in coinfection with R. viridis (37.72%) predominated. Also, 5.89% were singly infected with R. viridis and 8.64% were deprived of facultative endosymbiont. The α-diversity in aphid endosymbiont infection statuses varied according to the pea aphid biotype only (LM: F = 26.13, df = 2, p < 0.001): the Shannon index was the highest FIGURE 2 | Frequencies of the parasitoid species in the different pea aphid biotypes ("alfalfa," "clover," and "pea") at different parasitism period ("pre-parasitism," "maximal parasitism," and "post-parasitism").
in the alfalfa fields, intermediate in pea and the lowest in clover (Figure 3).
When we analyzed the relationship between the Shannon Index (H) of parasitoid communities and the Shannon Index (H) of endosymbiont infection statuses of aphids, we found a significantly positive covariation (Pearson correlation, r = 0.499, p = 0.004) in α-diversities (Figure 3).

β-Diversity in Parasitoid Communities and Endosymbiont Infection Statuses of Aphids
When the Bray-Curtis distances were calculated between all pairwise combinations of parasitoid communities, we found that they varied significantly in composition between A. pisum biotypes (Figure 5). PERMANOVA detected a significant effect of aphid biotype on the parasitoid community assemblages (F = 3.84, df = 2, p = 0.003) and this factor accounted for 21% of the variance in the data. Pairwise comparisons between biotypes showed that the structure of parasitoid communities emerging from pea crops differed from the two other crops. No temporal dynamics in β-diversity in parasitoid communities was found (PERMANOVA, period effect: F = 1.565, df = 2, p = 0.150).
The Bray-Curtis dissimilarities between endosymbiont infection statuses of aphids presented highly contrasted values and showed that they differed strongly between biotypes while being stable during the season within each biotype (Figure 5).
PERMANOVA confirmed this pattern since the biotype effect on infection statuses dissimilarity was highly significant (F = 79.08, df = 2, p = 0.001). This factor accounted for 83% of the variance in the data. No temporal dynamics in β-diversity of endosymbiont infection statuses of aphids was found (PERMANOVA, period effect: F = 740, df = 2, p = 0.596), confirming their stability throughout the cropping season.
A positive correlation between the two dissimilarity matrices was found (Mantel test: z M = 0.186, p = 0.002, see Supplementary  Figure 2). The dissimilarity between parasitoid communities was therefore correlated with the dissimilarity between endosymbiont infection statuses of aphids.

Parasitism Rates of Pea Aphids in Relation With Their Symbiotic Associations
In both clover and pea crops, the parasitism rate of all parasitoid species or of the dominant ones did not vary according to the prevalence of the most frequent symbiotic associations (Table 1 and Supplementary Figure 3). For the alfalfa biotype, while the prevalence of H. defensa in monoinfection did not influence the parasitism rates (Table 1 and Supplementary  Figure 3), the coinfection with H. defensa and F. symbiotica had a significantly negative effect on parasitism rates: when pea aphids from alfalfa fields presented high prevalence of coinfection with both protective endosymbionts, the parasitism rate of A. ervi and the overall parasitism rate declined although non-linearly (Table 1 and Figure 6).

Patterns of Temporal Variation in H. defensa Strain Diversity
A subsample of 72 A. pisum from different alfalfa fields and infected with H. defensa singly or co-infected with H. defensa and F. symbiotica was analyzed to characterize the strain variation of H. defensa. Among the pea aphid individuals, 36 aphids were selected before the peak of parasitism and 36 aphids after the peak. The APSE phage associated with H. defensa was consistently detected in all aphids. Seven genetically different strains of H. defensa were characterized ( Table 2). Two main haplotypes, representing 73% of the H. defensa strains, dominated aphid populations before and after parasitism peak. A greater diversity of strains was observed after the peak of parasitism (0.11 before and 0.19 after) and a similar pattern was observed for the aphid clonal diversity (0.58 before and 0.64 after). However, the frequencies of H. defensa strains did not differ significantly before and after the parasitism peak in alfalfa fields (χ 2 test: χ 2 = 6.65, df = 6, p = 0.645).

DISCUSSION
Since some microbial symbionts confer a protection against natural enemies to their hosts that can potentially alter food web interactions (Hafer and Vorburger, 2019;Mclean, 2019), our objectives were to test under natural conditions whether the parasitism rate of different A. pisum biotypes depended FIGURE 4 | The diversity and relative abundance of endosymbiont infection statuses of pea aphids in relation with aphid biotype ("alfalfa," "clover," and "pea") and parasitism period ("pre-parasitism," "maximal parasitism," and "post-parasitism"). For each column, the colored area corresponds to the percentage of aphid individuals harboring the corresponding symbiont species. When colored areas are present in several columns, it indicates a co-infection with multiple symbiont species. (Spi: Spiroplasma sp.; Reg: Regiella insecticola; Ham: Hamiltonella defensa; Rica: Rickettsiella viridis; Ser: Serratia symbiotica; Rick: Rickettsia sp.; Fuk: Fukatsuia symbiotica). on their endosymbiont communities and whether symbiont communities and parasitoid communities associated with pea aphids were related. We found that facultative endosymbiont communities were highly structured by biotype and stable in time while the parasitoid communities showed moderate differences between pea aphid biotypes and some change in structure over time. At the level of the pea aphid complex, we revealed a correlation between diversities (i.e., αand β-diversities) of endosymbiont infection statuses of aphids and parasitoid communities. Interestingly, we found a negative correlation between the prevalence of H. defensa and F. symbiotica in coinfection and the intensity of parasitoid pressure in the alfalfa biotype, confirming in field conditions the protective effect of this symbiotic combination.
The strong associations between endosymbiont communities and pea aphid biotypes are not novel and have been recurrently reported in various studies (Simon et al., 2003;Ferrari et al., 2004;Russell et al., 2013;Henry et al., 2015;Smith et al., 2015). Beyond the absence or the presence of each facultative bacterial endosymbiont in natural populations, each pea aphid biotype presented dominant symbiotic associations during the season: more than 80% of the aphids feeding on clover were singly infected with R. insecticola; in alfalfa, the infection with H. defensa alone or in coinfection with another facultative symbionts predominated the natural populations; and almost all aphids specialized on pea harbored S. symbiotica singly or in coinfection with R. viridis. Several hypotheses have been put forward for such differences in symbiotic associations and involve either ecological filters exerting selective pressure on aphid symbioses, symbiont-symbiont interactions or the effect of drift on symbiont associations (Mathé-Hubert et al., 2019). An essential step for linking facultative endosymbionts with parasitoid communities was to reliably assess the parasitism pressure exerted on aphid populations. Here, we used two proxies to reduce the bias of underestimating the parasitoid pressure by measuring only the mummification rate (Oliver et al., 2003). For the three biotypes, both parasitoid pressure proxies were highly Bold values refer to significant effect.
FIGURE 6 | Relationship between Hamiltonella defensa/Fukatsuia symbiotica coinfection and parasitism rate (assessed for the dominant parasitoid species, Aphidius ervi, and for all parasitoids) in pea aphids from alfalfa fields. ID refers to the alfalfa field. Solid line represents the predicted values from the fitted generalized additive mixed models.
positively correlated and showed the same temporal dynamics, indicating that the actual parasitoid pressure exerted on field populations of the pea aphid in different ecological situations was reliably assessed. Overall, the parasitism pressure over time presented a single peak of high parasitoids' activity in early July. This temporal dynamic was similar between alfalfa, clover and pea biotypes but varied quantitatively: the pressure exerted by parasitoids was the highest in clover, intermediate in pea and the lowest in alfalfa crops. One explanation for these quantitative differences could be that clover is typically grown in a more complex agricultural mosaic than the other crops, which would favor biological regulations. An alternative could be that clover is more attractive to parasitoids. Despite these temporal dynamics and inter-biotype variations in parasitism pressure, the relative abundance of the endosymbiont infection statuses of aphids changed very little over the season. These results contradict those obtained in a previous study done on both alfalfa and clover biotypes in the United States (Smith et al., 2015). In this earlier work, considerable seasonal shifts were indeed observed in the frequencies of endosymbionts, especially in Pennsylvania fields. This difference may result from the fact that this previous work considered each symbiont species individually while we analyzed them as communities. It could also be due to differences in population and community composition of endosymbionts and parasitoids, which exist between the native and introduced range of A. pisum. For example, we found six species of parasitoids in our survey while Smith et al. found only two. In addition, while there is evidence of some protection against parasitoids conferred by F. symbiotica in western Europe (Heyworth and Ferrari, 2015;Leclair et al., 2016), this does not appear to be the case in the United States (Doremus and Oliver, 2017), suggesting strain differences between countries. The temporal stability of symbiotic associations we observed suggests a good fidelity in vertical transmission of symbiont combinations in natural conditions, confirming previous results on field estimates of maternal transmission rate of pea aphid endosymbionts (Rock et al., 2018). This seasonal stability was also observed at the symbiont population level as we did not detect differences in the frequencies of H. defensa strains throughout the alfalfa growing season. Temporal fluctuations in symbiotic associations may, however, occur on a longer time scale (between years or beyond), although the same symbiotic associations have been found repeatedly in pea aphid populations in independent studies and in different years (Henry et al., 2013;Rock et al., 2018;Mathé-Hubert et al., 2019). Interestingly, we found a negative non-linear relationship between the frequency of co-infections with H. defensa and F. symbiotica and the parasitism rate (estimated for all parasitoid species or for A. ervi alone). Surprisingly, the single infection with H. defensa alone did not appear related to the rate of parasitism. However, such results are consistent with observations made in controlled conditions. Indeed, it has been shown that the protection conferred by H. defensa alone varies greatly according to the host and genotypes of H. defensa (Oliver and Higashi, 2019), whereas co-infections with F. symbiotica give the aphids a high or even total protection against parasitism (Guay et al., 2009;Leclair et al., 2016), but see Doremus and Oliver (2017). Overall, this endosymbiont protection may be responsible for the lower parasitism rate (between 30 and 50% reduction) observed in the alfalfa biotype compared to the two other biotypes. Negative correlations between protective endosymbionts in pea aphids and parasitoid rates have been reported in earlier field studies (Smith et al., 2015;Hrcek et al., 2016;Rothacher et al., 2016).
Our study and the previous ones thus confirm the effectiveness of symbiont-mediated protection in natural environments exploited by complex and diverse parasitoid communities.
One explanation given for the absence of H. defensa in clover and pea biotypes is because of a lower parasitoid pressure in natural populations of these biotypes (Oliver et al., 2008). However, our results showed that parasitoid wasps severely attacked these two pea aphid biotypes, with parasitism rates up to 70%. Given this parasitism pressure, an alternative protection against parasitoids may exist in these aphid populations; for instance, another facultative endosymbiont could confer an alternative protection to H. defensa. However, we did not observe a relationship between endosymbiont communities in the pea and clover biotypes and the rate of parasitism, although S. symbiotica (predominant in the pea biotype) and a strain of R. insecticola (the dominant symbiont in the clover biotype) have been reported to confer some parasitoid resistance in several laboratory studies (Oliver et al., 2003;Vorburger et al., 2010;Heyworth and Ferrari, 2015). Also, we showed in controlled conditions that the endosymbiont communities that dominated clover and pea biotypes in our field survey conferred very limited protection to A. ervi (Leclair et al., pers. obs.). Other ecological and non-ecological factors would better explain the prevalent symbiotic associations observed in these biotypes as discussed earlier (Mathé-Hubert et al., 2019). For example (Smith et al., 2015) suggested that the high prevalence of R. insecticola in the clover biotype was due to the strong incidence of fungal pathogen-induced mortality in clover fields, which would select this symbiont because of the fungal protection it confers (Scarborough et al., 2005). More field works are needed to assess the influence of environmental factors on symbiont composition of the various pea aphid biotypes.
Two earlier field studies showed that protective symbionts could influence the third trophic level by shaping the structure of parasitoid community attacking aphid populations (Hrcek et al., 2016;Rothacher et al., 2016). In addition, an experimental evolution study showed that parasitoid diversity could maintain diversity in protective symbionts (Hafer-Hahmann and Vorburger, 2020). Our study is in line with this body of work linking symbiont and parasitoid diversities. Indeed, we found that both αand β-diversities of symbionts and parasitoids were correlated, suggesting some interactions between these two communities through their aphid hosts or other environmental factors (i.e., local habitats).
Given the high prevalence of H. defensa in the alfalfa biotype and the negative correlation we found only between H. defensa-F. symbiotica co-infection and parasitoid pressure, one might have expected to find a strong effect of the symbiotic protection found in the alfalfa biotype on the community structure of parasitoids attacking the pea aphid. However, we showed that parasitoid communities varied little between the clover and alfalfa biotypes but more between these two and the pea biotype. A. ervi is often cited as the parasitoid species exerting the highest pressure on pea aphid populations (Kavallieratos et al., 2004;Smith et al., 2015). The alfalfa biotype was predominantly attacked by this parasitoid species despite the high prevalence of protective symbioses in this biotype. Previous works showed an evolutionary potential of parasitoids to counteract the symbiont protection conferred by H. defensa (Dion et al., 2011;Rouchet and Vorburger, 2014); some genotypes of A. ervi may thus have evolved virulence factors against H. defensa mediated protection in the wild (Dennis et al., 2020). Our results showed that the parasitoid community associated with the pea biotype strongly differed from the two others, in particular by the abundance of A. eadyi and A. avenae. Ferrari et al. (2004) have showed that H. defensa could provide resistance to A. eadyi in A. pisum populations. However, this endosymbiont was not present in the pea biotype in our survey and it is not known whether S. symbiotica or R. viridis, which, on the other hand, were both well represented in A. pisum from pea, confer resistance to A. eadyi. More work is needed to test whether the most prevalent symbiotic combination in each biotype confers an optimal protection to the corresponding parasitoid communities.
In conclusion, this study showed a temporal stability in symbiont populations and communities, in sharp contrast with a strong seasonality in parasitoid activities. This weak response of symbiont communities to parasitoid pressures could be explained by the limited costs of carrying protective symbionts on this timescale. However, further field surveys are needed to determine whether the composition in protective symbionts is maintained over longer period or is rather driven by an ecological-evolutionary dynamics resulting from selection for or against resistance, as suggested earlier (Smith et al., 2015) and recently demonstrated in a manipulated agricultural system (Ives et al., 2020). We also showed that the three pea aphid biotypes, despite their distinct endosymbiont composition, were exposed to similar range of parasitism pressures, suggesting that other protective alternatives than hosting H. defensa alone or with F. symbiotica, and involving symbiont or host mechanisms, may be used by the various biotypes of the pea aphid complex, which merits further work. Finally, we detected a link between communities of parasitoids and symbionts, suggesting interactions through shared resources or other environmental filters. The study of other communities of natural enemies of aphids (e.g., predators, pathogens) could reveal more such links and allow to better measure the importance of symbionts in food webs.

DATA AVAILABILITY STATEMENT
The datasets generated and analyzed in this article are available in Zenodo at https://doi.org/10.5281/zenodo.4548282.

AUTHOR CONTRIBUTIONS
ML, J-CS, and YO conceived and designed the field works and experiments and managed the field samplings. FM and ML performed both molecular studies and symbiont detection. CB performed species identification. YO analyzed the data. YO and J-CS contributed to the supervision of this study. All authors contributed critically to the drafts and gave final approval for publication.

ACKNOWLEDGMENTS
We are very grateful to all people who helped in aphid sampling, to Anne-Lise Boixel and Antoine Briand for their strong investment in field surveys and laboratory works during summer 2014, and to Bernard Chaubet for his help in parasitoid species determination. This work was supported by the French "Ministère de l'Enseignement Supérieur et de la Recherche, " the Plant Health and Environment division of INRAE and a grant from the French Research Agency (HMICMAC project 16-CE02-0014).

618331/full#supplementary-material
Supplementary Figure 1 | Location of the sampled fields. The color and the shape of the dots represent the legume species ("Alfalfa," "Clover," "Pea") on which pea aphids and parasitoids were collected.
Supplementary Figure 3 | Effect of the most prevalent symbiotic associations on the parasitism rate (dominant parasitoid species and all parasitoid together) in pea aphid individuals from the alfalfa, clover and pea biotypes. ID refers to the field crop sampled. For all these panels, no significant relationship was found.
Supplementary Table 1 | Sampling information: code, location and sampling frequency of the fields considered in our study. For a given legume crop, the field code assignment depends on the geographical location: we numbered from the westernmost to the easternmost field.
Supplementary Table 2 | Sampling information: dates, field identity (i.e. field), number of aphid nymphs used to assess symbiotic composition (i.e. nymph), number of adult aphids kept to measure parasitism rate (i.e. adult) and its value for each sampling (i.e. para) (when "-", no sampling).