Spatial and Temporal Variation in the Aphid–Parasitoid Interaction under Different Climates

in the Aphid–Parasitoid Interaction Abstract: Global warming will increase pest insect population sizes and diminish the effectiveness of biological control. This biological control failure scenario appears to be of particular concern for areas with a signiﬁcant increase in maximum temperatures, such as the increase experienced in the Central Valley of Chile over the last 40 years. We assessed the impact of different climatic zones and maximum temperatures along the coast and the Chilean Central Valley on the grain aphid ( Sitobion avenae ) density, parasitism rate, and facultative endosymbionts in wheat ﬁelds during the growing season in the springs of 2017 and 2018. A signiﬁcant effect on aphid density due to zones and maximum temperatures was detected; however, this depended on the zone and year analyzed. Changes between zones and seasons were observed for parasitism rates, while maximum temperatures only signiﬁcantly affected the parasitism rate in 2017. The main parasitoid wasp found was Aphidius ervi in both zones and seasons. Regiella insecticola infected 95% of the samples in both zones, although it does not seem to have a protective role at the ﬁeld level. Our ﬁndings suggest that, at present, global warming does not signiﬁcantly affect the grain aphid outbreaks and their biological control in Chile. However, this study points out the importance of pre-emptive monitoring to detect aphids and the synchrony loss of their parasitoid wasps.


Introduction
Climate change (e.g., increased environmental temperature, increased atmospheric CO 2 , unstable climates, and altered frequency/intensity of extreme weather events) is the most serious concern for agriculture [1,2]. Global warming, characterized by increased and extreme temperatures that will become more frequent, affects crop yields, mainly through crop pest biology and distribution changes, especially in invasive pest species [3][4][5].
Aphids (Hemiptera: Aphididae) are a highly invasive pest species due to their broad phenotypic plasticity, which includes (a) their reproduction mode (see [6]), (b) the systematic development of insecticide resistance, (c) the development of defenses against plant chemistry and natural enemies, and (d) their symbiosis with obligate and facultative bacteria, which confer aphids with several relative advantages, including protection against natural enemies and heat tolerance [7,8]. A worldwide pest of cereals, the English grain aphid, Sitobion avenae (Fabricius) [9,10], originated from Europe and the Mediterranean zone and was introduced to Chile in the late 1960s, causing significant losses to cereal production [11]. Mild winters in Chile favored asexual reproduction all year round, allowing aphid populations to distribute rapidly and widely across most of the cereal production area [12]. Several Aphidiinae endoparasitoid species were introduced to Chile from different geographic origins during the late 1970s to cope with this pest as part of a governmental biocontrol program [13]. Hence, a parasitoid assemblage (i.e., the set of parasitoid species attacking aphids) has provided an essential contribution to the suppression of high population densities of cereal aphids in Chile to date, with minimal input from insecticides [14]. Indeed, Chile is considered one of the best examples of successful aphid biological control programs to combat the English grain aphid [15].
Global warming, however, threatens the dynamics, synchronization, and temporal and spatial structures of the interactions among pests and their natural enemies [16]. Global warming can potentially increase the frequency of bacterial endosymbionts that provide aphids with abiotic protective effects, such as heat tolerance, even further [17][18][19][20][21]. If their transmission and geographic distribution increase at higher temperatures, the success of the biological control program may be threatened [21][22][23][24].
Of particular concern are the maximum daily temperatures and frequency of warm events. These variables have the most critical effect on the life history traits of insect pests and biological control agents, making pest-derived problems likely to become more unpredictable [25][26][27]. The troubling fact is that Chile's central valleys have shown one of the most striking increases in maximum temperatures in recent decades, which threatens to continue [28,29].
The monitoring of the temporal dynamics of interacting species in the agroecosystems is crucial to anticipate climate change effects on aphid outbreaks and biocontrol populations, especially those prone to global warming, highlighting the urgent need for these types of studies [16,30]. Hence, this study aimed to determine whether spatial and temporal variations in maximum temperatures modify the density of both S. avenae aphids and their parasitism rate in the field. We assessed the aphid density, parasitism rates, and frequency of facultative endosymbionts over two consecutive years throughout the wheat growing season, comparing the situation among field crops located in two areas-one with temperatures regulated by the Pacific Rim and the other with significant temperature increases in recent decades. We discuss the results in terms of management practices that ensure the continuity of the biological control of S. avenae.

Sampling Sites and Climates
Surveys were conducted in commercial insecticide-free spring wheat fields in two zones with contrasting climatic conditions (Figure 1) (Table S1, Supplementary Materials). Zone 1 is located close to the Pacific Rim and displays a coastal temperate Mediterranean climate, with rainfed wheat crops. Zone 2 is in the middle of the Chilean Central Valley, where a warmer Mediterranean temperate climate dominates, and wheat crops need periodic irrigation (four irrigation periods per season, occurring from September to December). The distance between fields in each zone ranged from 1.5 km to 19 km. Because of crop rotation, some wheat fields were different across the two years.
Meteorological data were collected from the national agrometeorological network [31]. We selected wheat fields located close to a meteorological station to achieve more accurate environmental temperature measurements. Hence, one meteorological station was used per zone (two in total). Figure S1 shows the temperatures recorded on each sampling date. We computed the average annual temperatures from the last three years of data (2015-2018) to develop a clearer picture of the climatic differences between the studied zones. Thus, Zone 1 s temperatures ranged between 17.5 • C (max) and 7.8 • C (min) with 742 mm of rainfall, while Zone 2 s temperatures ranged between 21 • C (max) and 7.3 • C (min) with 668 mm of rainfall. Moreover, Zone 2 showed an alarming increase in its maximum temperature in recent decades [29,32].  (Table S1).
Meteorological data were collected from the national agrometeorological net [31]. We selected wheat fields located close to a meteorological station to achieve accurate environmental temperature measurements. Hence, one meteorological s was used per zone (two in total). Figure S1 shows the temperatures recorded on each pling date. We computed the average annual temperatures from the last three ye data (2015-2018) to develop a clearer picture of the climatic differences between the ied zones. Thus, Zone 1′s temperatures ranged between 17.5 °C (max) and 7.8 °C with 742 mm of rainfall, while Zone 2′s temperatures ranged between 21 °C (max) an °C (min) with 668 mm of rainfall. Moreover, Zone 2 showed an alarming increase maximum temperature in recent decades [29,32].
Because both zones had similar average temperatures, we determined the num >30 °C events and the average maximum temperatures (2015-2018) during the a spring (September-December), which is the growing period for most wheat cultiv Chile, as indicators. Zone 1 showed no days with temperatures >30 °C and an av maximum temperature of 17.5 °C, while Zone 2 showed 37 days with temperature °C and an average maximum temperature of 22.6 °C.

Aphids and Parasitoids Samplings
Three wheat fields were sampled in each zone. To standardize the sampling e we collected aphids and parasitized aphids (referred to as "mummies") from a prede area of 1 ha per field, avoiding borders (Table S1). Live aphids and mummies wer lected every 14 days from the wheat tillering stage to the dough-ripening stage, from tember to December (austral spring), for two consecutive years (2017 and 2018). We pled the total population of aphids and mummies observed on 20 wheat tillers in randomly selected sampling points in each field, separated by at least 3 m from the border and 10 m from one another [33,34]. The leaves and ears of the tillers infected aphids were clipped off and carefully transferred to modified Petri dishes with venti  (Table S1).
Because both zones had similar average temperatures, we determined the number of >30 • C events and the average maximum temperatures (2015-2018) during the austral spring (September-December), which is the growing period for most wheat cultivars in Chile, as indicators. Zone 1 showed no days with temperatures >30 • C and an average maximum temperature of 17.5 • C, while Zone 2 showed 37 days with temperatures >30 • C and an average maximum temperature of 22.6 • C.

Aphids and Parasitoids Samplings
Three wheat fields were sampled in each zone. To standardize the sampling effort, we collected aphids and parasitized aphids (referred to as "mummies") from a predefined area of 1 ha per field, avoiding borders (Table S1). Live aphids and mummies were collected every 14 days from the wheat tillering stage to the dough-ripening stage, from September to December (austral spring), for two consecutive years (2017 and 2018). We sampled the total population of aphids and mummies observed on 20 wheat tillers in five randomly selected sampling points in each field, separated by at least 3 m from the field border and 10 m from one another [33,34]. The leaves and ears of the tillers infected with aphids were clipped off and carefully transferred to modified Petri dishes with ventilation and sealed with Parafilm ® .
All fields were sampled on the same day and repeated every two weeks to complete six sampling dates (T1-T6). Once in the laboratory, aphids and mummies were counted and identified following taxonomic keys [9]. Mummies were kept in Petri dishes until the emergence of wasps, and the parasitoids were identified using taxonomic keys [11]. Live aphids were kept in cages containing wheat seedlings (Triticum aestivum cultivar Pantera) for 15 days and checked for mummification every two days. New S. avenae mummies were removed from the seedlings and kept in Petri dishes for counting and identification.
We collected an additional sample of 10-30 adult S. avenae aphids during each sampling date (field, zone, and year). We collected aphids through a linear transect to limit the chances of sampling individuals that belonged to the same parthenogenetic colony (i.e., clonal lineages) by taking one single wingless individual every 10 m. All aphids from those samples were stored in Eppendorf tubes with 95% ethanol until bacterial endosymbiont screening (700 aphids in total) (Table S2).

Aphid Bacterial Endosymbionts
The total DNA from each stored aphid was extracted by the salting-out method [35]. We screened for the presence of the most frequent endosymbionts harbored by aphids (Hamiltonella defensa, Regiella insecticola, Serratia symbiotica, Fukatsuia symbiotica, Rickettsia sp., Ricketsiella sp., and Spiroplasma sp.) using a polymerase chain reaction based protocol (PCR) described previously [36,37]. The obligate endosymbiont, Buchnera aphidicola, was used as a positive control.

Data Analysis
All data were analyzed in the R software version 3.6.1. [38]. Differences in maximum temperatures were log-transformed and analyzed with a linear model.
To test the effect of every zone, year, and sampling date on S. avenae, parasitism rate, and facultative endosymbiont proportion, we used generalized linear mixed models (GLMMs).
The density of the English grain aphid S. avenae was analyzed as the number of individuals per 100 tillers in each field, with a Poisson distribution error. Due to the low number of emerged parasitoids in some fields and sampling dates, we did not perform statistical analysis for the composition of parasitoid species; therefore, data are only descriptive and were calculated as the proportion of each identified parasitoid species to the total of emerged parasitoids for every zone, year, and sampling date. The parasitism rate was calculated only for those fields with a frequency of S. avenae >6. Hence, the parasitism rate was calculated as the proportion of parasitized S. avenae individuals, i.e., the number of mummies found in 100 tillers/number of aphids plus the number of mummies found in 100 tillers [33] for every zone, year, and sampling date, and analyzed with a binomial distribution error. Changes in the density of S. avenae during the season were evaluated for every zone and year.
Because R. insecticola was the predominant facultative endosymbiont (except for a single aphid individual bearing a co-infection with R. insecticola and H. defensa), we calculated the proportion of infected aphids as the number of aphids carrying R. insecticola divided by the total number of individuals collected per year, zone, and sampling date (Table S2), analyzed with a binomial distribution error.
We used the year, zone, and sampling date (T1-T6) as fixed parameters and the field in each zone as a random factor. The effect of each factor was evaluated by type II Wald chisquared tests [39]. Comparisons of each fixed factor for S. avenae, parasitism rate, and facultative endosymbiont proportion were performed with the R package multcomp [40]. Because maximum temperatures differed between zones and years, the effect of maximum temperature on S. avenae, parasitism rate, and facultative endosymbiont proportion was analyzed separately for every zone and year using a GLMM with the maximum temperature as a fixed factor and the field as a random factor.

Registered Temperatures in the Field
Our results show that Zone 2 registered higher maximum temperatures than Zone 1 in both seasons across the sampling dates (2017: F (F-statistic) = 93.2, df = 1, p ≤ 0.001; 2018: F = 68.5, df = 1, p ≤ 0.001) ( Figure S1). In 2017, Zone 1 displayed a maximum temperature range of 11.8-21.2 • C, with an average maximum temperature of 16.1 • C, while, in 2018, it displayed a range of 11.7-21.5 • C, with an average of 17.2 • C. No extreme maximum temperature events were recorded for either of the studied seasons (>30 • C). Zone 2 displayed a maximum temperature range of 11.6-31.1 • C in 2017 (with an average maximum temperature of 21.3 • C) and one single extreme high-temperature event (between T5 and T6), and 13.2-32.7 • C in 2018, with an average maximum temperature of 22 • C and two extreme high-temperature events (between T3-T4 and T5-T6). Data for the daily accumulated rainfall are provided in Figure S2.

Sitobion avenae Density
For the 2017 campaign, we collected 1305 aphids and mummies, from which 68.3% were identified as S. avenae. The S. avenae sample contained 192 mummies. For the 2018 campaign, we sampled 1128 aphids, 84.4% of which belonged to S. avenae, with 123 mummies (Table S2). We found a maximum of 109 S. avenae aphids per 100 wheat tillers (1.09 aphids per tiller on average).
When the effects of the zone, year, and sampling date on S. avenae density were tested, the results showed that S. avenae did not change significantly between the two years. However, the zone, the sampling date, and the interaction between the three factors had a significant impact ( Table 1). The differences between zones were mainly due to the fact that there were significantly lower aphid densities in Zone 1 than Zone 2 (Table S2). In contrast, differences in zones by year were due to a decrease in the aphid population in Zone 1 from 2017 (n = 495) to 2018 (n = 343) (Table S2). Chi-squared (χ 2 ), degrees of freedom (df), and p-value (p) for each effect (and their interactions) in the type II Wald chi-squared tests performed for the Generalized linear mixed models.
Changes in the density of S. avenae during the season were analyzed separately for every zone and year. In 2017, a higher number of aphids was observed in Zone 1 at the beginning of the season (χ 2 = 20.94, df = 5, p = 0.0008) (T1; Figure 2A). However, this effect was not observed in Zone 2 (χ 2 = 6.08, df = 5, p = 0.30) ( Figure 2B). In 2018, the density of S. avenae also changed throughout the season in Zone 1 (χ 2 = 17.03, df = 5, p = 0.004), where a lower number of aphids was observed at T3 compared to T6 ( Figure 2C). Similarly, the density of the S. avenae also changed throughout the season in Zone 2 (χ 2 = 22.8, df = 5, p ≤ 0.01), with a marked decrease at T6 ( Figure 2D).

Parasitoid Composition
Seven parasitoid species belonging to the Braconidae family were found. Six spec belonged to the Aphidius genus (Aphidius ervi, A. uzbekistanicus, A. rhopalosiphi, A.

Parasitism Rate
The zone had no effect on the parasitism rate; however, the year parameter had a small but significant effect. The sampling date and the interactions among the three factors were also significant ( Table 1). A significant increase in the parasitism rate was observed in Zone 1 from 2017 to 2018 (χ 2 = 1948.2, df = 1, p ≤ 0.001; 17.19% and 20.65%, respectively), while the opposite occurred in Zone 2 (χ 2 = 6.2623, df = 1, p = 0.012) (20.78% and 13.05%, respectively, for 2017 and 2018) (Table S2).
Changes in the parasitism rates at different sampling dates were analyzed separately for every zone and season. The sampling date showed a significant effect on the parasitism rate in 2017 for Zone 1 (χ 2 = 43.57, df = 5, p ≤ 0.001), with most differences observed when T1 (which displays a lower parasitism rate) was compared with T4 and T6 (Figure 2A). In contrast, the parasitism rate did not change throughout the season in Zone 2 (χ 2 = 7.91, df = 5, p = 0.16) ( Figure 2B). In 2018, very low S. avenae densities were found in most of the fields belonging to Zone 1 (Table S2), making it impossible to collect enough data to calculate the effect of the sampling date on the parasitism rate in this season; therefore, we did not include comparisons for this season, but the data are still shown in Figure 2C.
No variation in the parasitism rate of S. avenae was detected throughout 2018 in Zone 2 (χ 2 = 6.74, df = 5, p = 0.240). Because low aphid densities were detected at T3 and T6 in most of the fields in this zone, we did not calculate parasitism rates for these sampling dates ( Figure 2D), and they were not included in the multiple comparisons.

Parasitism Rate
The zone had no effect on the parasitism rate; however, the year parameter had a small but significant effect. The sampling date and the interactions among the three factors were also significant ( Table 1). A significant increase in the parasitism rate was observed in Zone 1 from 2017 to 2018 (χ 2 = 1948.2, df = 1, p ≤ 0.001; 17.19% and 20.65%, respectively), while the opposite occurred in Zone 2 (χ 2 = 6.2623, df = 1, p = 0.012) (20.78% and 13.05%, respectively, for 2017 and 2018) (Table S2).
Changes in the parasitism rates at different sampling dates were analyzed separately for every zone and season. The sampling date showed a significant effect on the parasitism rate in 2017 for Zone 1 (χ 2 = 43.57, df = 5, p ≤ 0.001), with most differences observed when T1 (which displays a lower parasitism rate) was compared with T4 and T6 (Figure 2A). In contrast, the parasitism rate did not change throughout the season in Zone 2 (χ 2 = 7.91, df = 5, p = 0.16) ( Figure 2B). In 2018, very low S. avenae densities were found in most of the fields belonging to Zone 1 (Table S2), making it impossible to collect enough data to calculate the effect of the sampling date on the parasitism rate in this season; therefore, we did not include comparisons for this season, but the data are still shown in Figure 2C.
No variation in the parasitism rate of S. avenae was detected throughout 2018 in Zone 2 (χ 2 = 6.74, df = 5, p = 0.240). Because low aphid densities were detected at T3 and T6 in most of the fields in this zone, we did not calculate parasitism rates for these sampling dates ( Figure 2D), and they were not included in the multiple comparisons.

Composition of Facultative Endosymbionts
Out of the seven facultative bacterial endosymbiont species screened, only cola and H. defensa were detected in both years and zones (Table S2). The bacteri ticola was carried by 95% of the sampled aphids, mostly as single infections, wh the aphids showed no infection. Co-infections with R. insecticola and H. defens were observed in one individual, which was excluded from the analysis to avo tions. When the frequency of infections with R. insecticola was analyzed, no eff year, zone, or sampling date were detected ( Table 1)

Discussion
Climate can shape species' distribution and their interactions, and temperat of the most critical climatic factors. Changes in environmental temperatures may impact the geographic distribution of pests and their biological control agents [ The effect of high temperatures has been broadly studied, mainly under la controlled conditions. The potential impact of daily maximum temperatures, has been poorly studied at the field level. Because the daily maximum temper increased rapidly in recent years and has been shown to have a more substantia pest insects and their interacting species [42], it is a more realistic approach for e global warming's effects on biological control. Our study considered two zones ferences in their maximum temperatures, and we assessed the impact of these the temporal densities of S. avenae and their parasitoids. Furthermore, we consi frequency of infections with facultative endosymbionts in S. avenae, which are ex have a role in the adaptive potential of aphids to face global warming and the tions with parasitoids [42].

Composition of Facultative Endosymbionts
Out of the seven facultative bacterial endosymbiont species screened, only R. insecticola and H. defensa were detected in both years and zones (Table S2). The bacteria R. insecticola was carried by 95% of the sampled aphids, mostly as single infections, while 4.9% of the aphids showed no infection. Co-infections with R. insecticola and H. defensa bacteria were observed in one individual, which was excluded from the analysis to avoid distortions. When the frequency of infections with R. insecticola was analyzed, no effects of the year, zone, or sampling date were detected ( Table 1). The maximum temperatures had no effect on R. insecticola proportions in any zone and year (Zone 1, 2017: χ 2 = 2.1, df = 1, p = 0.145; Zone 2, 2017: χ 2 = 0.4, df = 1, p = 0.517; Zone 1, 2018: χ 2 = 0.02, df = 1, p = 0.884 and Zone 2, 2018: χ 2 = 0.003, df = 1, p = 0.952).

Discussion
Climate can shape species' distribution and their interactions, and temperature is one of the most critical climatic factors. Changes in environmental temperatures may therefore impact the geographic distribution of pests and their biological control agents [41].
The effect of high temperatures has been broadly studied, mainly under laboratorycontrolled conditions. The potential impact of daily maximum temperatures, however, has been poorly studied at the field level. Because the daily maximum temperature has increased rapidly in recent years and has been shown to have a more substantial effect on pest insects and their interacting species [42], it is a more realistic approach for evaluating global warming's effects on biological control. Our study considered two zones with differences in their maximum temperatures, and we assessed the impact of these factors on the temporal densities of S. avenae and their parasitoids. Furthermore, we considered the frequency of infections with facultative endosymbionts in S. avenae, which are expected to have a role in the adaptive potential of aphids to face global warming and their interactions with parasitoids [42].

Factors Shaping Spatial and Temporal Sitobion avenae Populations
Our results show a clear difference in S. avenae densities between zones, with a higher aphid density in the warmer zone (Zone 2). Even if temperatures above 30 • C reduce the fecundity and extend the physiological development period of pests, reducing their population growth [43], an increase in temperature, as long as it does not exceed the thermal threshold of a given species, has been shown to have positive effects on aphid populations [44]. In S. avenae, significantly higher population growth is reported at 25 • C than at 15 • C [45]. Therefore, the temperatures recorded in Zone 2 seem to have a positive impact on S. avenae population numbers, which explains the higher density of S. avenae in that zone. It is important to note that Zone 2 also shows some extremely high-temperature events (>30 • C). However, they do not seem to have negatively affected S. avenae ( Figure S1) because no decrease in S. avenae density was observed after these events. The frequency of extremely high temperature events can negatively impact aphid biology because aphids may not have enough time and resources to recover from the potential heat damage produced by these events [46]. For example, field experiments in China showed a decrease in S. avenae densitywhen the frequency of extreme temperatures (>30 • C) was artificially increased by 60% compared with natural conditions [27]. Intriguingly, in our study, maximum temperatures during the season only positively affect aphid density in Zone 1 for 2018 (from T4 to T6). As we previously described, Zone 1 is located close to the Pacific Ocean, within a coastal, temperate Mediterranean climate with strong oceanic regulation. The increases in maximum temperatures observed in 2018 did not exceed the thermal threshold of S. avenae [45], instead, having a positive effect on S. avenae density. No effect of maximum temperatures was detected in Zone 2, probably because the maximum temperatures recorded throughout the seasons under study were not extremely high but close to the optimal temperature that would allow for an S. avenae population increase.
Although our study was focused on the climate particularities of each zone, mainly temperatures, it is important to keep in mind other factors, which were not studied, that could also influence temporal variations in S. avenae density(i.e., host plant phenology, the effects of biological control agents, landscape, rainfall, etc.). Temporal changes in S. avenae densityamong sampling dates were detected for both zones and seasons. Our study considered zones with different climates and temperatures that could affect aphid performance and their host plants' phenology [47]. In cereals, the highest peak of aphid density is usually observed at the milk-ripening stage [33], but the timing of this crop growth stage may vary according to the climate. Moreover, the synergy, additive, or antagonistic interaction of natural enemies, including parasitoids, predators, and pathogens (i.e., viruses, bacteria, and fungi) can decrease aphid density in fields and account for temporal variations in aphids throughout the season [48][49][50]. The complexity of the landscape could also impact aphid density because complex landscapes contain various grasses that could act as migration sources, thereby transferring aphids to wheat fields [51]. Zone 2 is characterized by a higher area of arable land than Zone 1; however, we did not find a lower aphid density in this zone, as expected. Hence, the landscape did not appear to be an important factor in explaining the S. avenae densities found in our results.
On the other hand, climate change and global warming favor drought in many areas worldwide [52], which may indirectly impact aphid density [53][54][55]. For instance, water deficits can increase the concentration of nitrogen and amino acid contents in the host plants of aphids, thus benefiting their growth and reproduction [56,57]. Moreover, water deficits can also decrease the parasitism rate by altering the parasitoid preference, thus impacting aphids' biological control [58]. Since 2010, central Chile has faced a dramatic sequence of dry years, with deficits in rainfall between 25% and 45%, referred to as the "Mega Drought," characterized by a few events of extreme precipitation followed by long periods with no rain [59]. Because Zone 2 has higher temperatures than Zone 1, this must mean that it bears the consequences of the drought of the last decade more significantly, thereby increasing the water deficit on plants, which could also account for the higher density of S. avenae in this zone.

Parasitoid Assemblage
Most parasitoid wasps collected in wheat fields belonged to the sub-family Aphidiinae, in which Aphidius ervi represented 82% of the whole sample. This observation agrees with previous reports showing that A. ervi and A. uzbekistanikus are the most common parasitoids attacking cereal aphids in Chile [50]. Interestingly, the parasitoid assemblage changed as the season progressed. While the most frequent parasitoid in both zones and seasons was A. ervi early in the season, other parasitoid species gained relevance later; nevertheless, these results must be interpreted with caution because very low numbers of parasitoids emerged at some sampling dates. Changes in population dynamics throughout the season seem to be a feature of the parasitoid assemblages that parasitize cereal aphids, as previously observed in various European countries [60][61][62]. The seasonal activities of some parasitoid species seem to match their thermal tolerances. Some studies [63,64] have reported that the thresholds for A. ervi from egg to mummy were 2.2 • C and 6.6 • C for mummy to adult development, respectively, while A. rhopalosiphi were 4.5 • C and 7.2 • C, and those for P. volucre were 3.8 • C and 5.5 • C, respectively. Early parasitization is critical in preventing aphid outbreaks [60,65]. A parasitoid able to develop at low temperatures could mean an earlier appearance in the field, which is essential for parasitizing aphids, as they start to multiply early. This feature does not appear to be affected by the climate conditions at present and may explain the early presence of A. ervi in wheat fields found in our study and the biocontrol of S. avenae in Chile. Nevertheless, permanent temporal and spatial studies similar to ours are needed to monitor the biological control status for this pest and others in the coming years, as we discuss below.

Spatial and Temporal Changes in Parasitism Rate
The demographic balance between aphids and parasitoid wasps is critical for the efficacy of biological control. Nevertheless, this balance can suffer from increased environmental temperatures due to phenological desynchronization between aphids and parasitoids, which disrupts their trophic relationships and ultimately contributes to aphid outbreaks [5,66].
Our results show that the parasitism rate varied between seasons and zones. Zone 1 showed a higher parasitism rate in 2018 than in 2017, while the opposite trend was observed in Zone 2. However, these results should be interpreted with caution since we had to remove data with a low number of aphids to avoid overestimating parasitism rates (Table S2). It has been reported that A. ervi shows a better parasitism rate at 25 • C than at 15 • C or 20 • C and the worst performance at 30 • C [67]. The rise in the maximum temperature from 2017 to 2018 in Zone 1 seemed to positively impact the parasitism rate because the colder zone (Zone 1) became a bit warmer ( Figure S1). But this produced a contrasting effect in the warmer zone (Zone 2), where temperatures are continuing to rise to levels at which the performance of parasitoid wasps significantly decays [5,66,67]. It is well known that species at higher trophic levels, such as parasitoids, are highly susceptible to temperatures close to 30 • C [42,68]. Similar to those recorded for Zone 2 in 2018, high temperatures can decrease the parasitoid attack rate [69] and alter life-history traits, such as the developmental time and lifespan of parasitoids [41]. Hence, during the last decade in central Chile, global warming-related phenomena could be following a sustained and troubling trend [29], which explains the differences in the parasitism rates observed between zones and seasons in this study and alerts us to the outcomes of biological control programs in upcoming years. More monitoring is needed to test this hypothesis.
Intriguingly, the maximum temperatures recorded throughout the season only influenced the parasitism rate in 2017. Assessing the effect of temperature on the parasitism rate at the field level is challenging because a single species does not attack aphids, but an assemblage of parasitoid wasps. Several parasitoid species with different thermal tolerances can attack S. avenae aphids; therefore, the outcome of parasitism depends on each parasitoid species' thermal tolerances. A more random parasitism rate was found in 2018 according to maximum temperatures and sampling dates (which were higher than those in 2017); this could mean that temperatures may affect each parasitoid species differently. For instance, the parasitoid wasp A. avenae is more tolerant of high temperatures, whereas A. rhopalosiphi is more susceptible [63], while A. ervi seems to tolerate low temperatures, but not temperatures over 30 • C [67,70]. Because A. ervi is the primary parasitoid attacking S. avenae in Chile, the predicted global warming could result in pernicious consequences for the biological control of S. avenae and wheat yields.
It is important to consider other factors that could impact parasitism rates, such as secondary parasitoids and predator attacks. Secondary parasitoids can cause a mortality rate of over 50% in parasitoid larvae [71], while predators can alter the parasitism rate by consuming parasitized aphids and adult wasps [72]. We only observed the appearance of secondary parasitoids at the last sampling date; hence, they probably did not impact our results significantly. On the other hand, predators were rarely observed during the samplings; however, their effect should not be entirely excluded.

Facultative Endosymbionts
Global warming can also affect the insect microbiome because increased temperatures have been shown to alter the composition of bacterial endosymbionts, which mediate aphid-parasitoid interactions [68,73].
All aphids carry an obligated bacterial endosymbiont, B. aphidicola, which provides aphids with essential amino acids and vitamins that they cannot obtain from plant sap [74]. Moreover, aphids can harbor various facultative endosymbionts that are non-essential but can confer relative advantages to aphids [75]. Facultative endosymbionts can increase thermal tolerance in aphids [18], broaden the host plant range [76], and offer protection against natural enemies (e.g., parasitoids and predators) [77][78][79]. Facultative endosymbionts are considered a significant phenotypic variation source, mainly in introduced aphid pest species, which are frequently characterized by asexual populations and reduced genetic diversities [7]. Hence, aphids can rapidly evolve in new environments and face global warming [42], which also has important implications for their interactions with their natural enemies [77,78,[80][81][82]. In warmer and dryer areas, Chilean populations of S. avenae are typically infected with R. insecticola, as the most frequent endosymbiont in central Chile [14,37]. Our results show a high frequency (95%) of R. insecticola in all fields, zones, and on all dates sampled, supporting the previously reported high frequencies of this bacteria, regardless of the range of maximum temperatures recorded in the two studied zones. It is likely that, because of the differences in both zones' temperatures, these temperatures are not sufficient to affect the facultative endosymbionts of S. avenae.
Increased temperatures can affect the frequency of infections with facultative endosymbionts. For instance, the rate of vertical transmission of R. insecticola indicates failure at >28 • C. However, when high temperatures are maintained for several generations, the transmission rate of R. insecticola is increased by up to 100%, showing that R. insecticola can improve aphid responses to heat [21]. Facultative endosymbionts providing environmental adaptations, such as tolerance to increased temperatures, could improve their prevalence quickly, which may increase and/or modify the distribution of aphids (e.g., in warmer climates) [18,83]. Whether the R. insecticola strain(s) found in our study provide aphids with heat tolerance should be tested under laboratory conditions. However, the high prevalence of strain(s) of R. insecticola in the field during the wheat growing season is an exciting observation, suggesting that there is no cost in relation to the fitness of S. avenae and no evidence of a strong selection against non-infected aphids.

Conclusions
Increases in the environmental temperatures can significantly improve aphid performance, positively impacting aphid population numbers [21]. Adverse effects on parasitoid wasps' performance can disable biological control [68]. In classical biological control scenarios, when both the pest and its natural enemy are introduced species, they usually have reduced genetic diversity and limited evolvability due to the actions of evolutionary factors during their introduction in the form of founding effects, selection (e.g., insecticide applications, natural enemies, climate), and genetic drift (e.g., reduced population sizes) [12,84,85]. Therefore, environmental stability appears to be an essential factor that maintains the balance of population sizes threatened by global warming. It can alter the herbivores in natural and managed systems and change the interactions between aphids and their natural enemies [86].
Nevertheless, at present, the biological control of S. avenae by parasitoid wasps does not appear to be significantly affected in Chile since the densities of S. avenae are still below the threshold of economic damage (five aphids per tiller [87]). Furthermore, facultative endosymbionts such as R. insecticola do not appear to have any protective effect against parasitism at the field level [14]. However, noteworthily, the warmer year in our study produced lower parasitism rates in the warmer zone, which could already be an indicator of the effects of global warming. Our research presents the current situation regarding the biological control of the English grain aphid in Chile at the field level, which may be useful for future comparisons. Finally, we suggest that periodic monitoring should be carried out to assess the densities of aphids and their facultative endosymbionts, as well as parasitoids and the inclusion of other biotic (e.g., hyperparasitoid and predator population dynamics) and abiotic (e.g., landscape composition, rainfall, etc.) factors that shape aphid-parasitoid interactions.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agriculture11040344/s1, Figure S1: Daily temperatures recorded for every zone and year throughout the season, Figure S2: Daily accumulated rainfall for every zone and year throughout the season, Table S1: Geographical coordinates and area of the wheat fields sampled in both studied zones and years, Table S2

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.