The epidemiology of BYDV is complex, involving both cultivated and non-cultivated grasses and multiple aphid vector species. In North American cereal agroecosystems, the cycle of disease begins when aphids colonize newly emerged winter wheat seedlings in the fall. Both viruses and vectors overwinter in wheat fields, but disease symptoms do not manifest until spring. Vector populations generally increase throughout the spring and then migrate to alternate grass hosts following wheat harvest. Summer cereals (e.g., corn and sorghum) and native perennial grass species are therefore an important “green bridge” maintaining the stability of the BYDV–wheat pathosystem (Brown, 1984; Clement et al., 1986; Ingwell and Bosque-Pérez, 2015). A suite of at least 25 aphid species transmit BYD Luteoviruses (Miller and Rasochová, 1997), with the predominant North American vectors being R. padi, S. avenae, and Schizaphis graminum (Gray et al., 1998; Hesler et al., 2005; Power et al., 2011). Transmission efficiency among vector species varies between BYDV strains and can be highly vector specific (Lucio-Zavaleta et al., 2001; Papura et al., 2002; Power et al., 2011). In general, transmission of Luteoviruses by aphids from plant to plant occurs via a persistent circulative process, whereby virus particles ingested during phloem feeding circulate through the aphid body prior to inoculation of a new host plant (Gray and Gildow, 2003; Ali et al., 2014).

We used a subset of English grain aphid (S. avenae) samples originally collected by Enders et al. (2018) across Kansas from winter wheat fields in 2015. All samples were originally stored at −80C directly after field collection and species identification. We selected samples collected from 6 wheat fields where BYDV had previously been detected and there were at least 10 individual aphids remaining in storage. Groups of 5 adult S. avenae were tested by Enders et al. (2018), thus we did not have an initial estimate of virus incidence (i.e., % viruliferous aphids) in these vector populations to aide in sample selection, only whether BYDV was detected. We extracted DNA and RNA from individual aphids using the AllPrep DNA/RNA Mini Kit (Qiagen) following manufacturer protocols, except that final elution volumes were 20 μl for DNA and 10ul for RNA. This dual DNA/RNA extraction allows for downstream DNA-based microbiome profiling (e.g., qPCR or metabarcoding) and screening for BYDV (RNA virus).

A second data set of Bird cherry-oat aphids (R. padi) was collected from winter wheat fields from 5 locations at Purdue Agricultural Centers across Indiana in 2018 and 2019. Weekly collections of aphids were performed using a standardized number of sweeps per field (100 sweeps/field) and timed with aphid population peaks (May–June), similar to methods used in Enders et al. (2018). Aphid samples were sorted, separated and identified to species, and then stored at −80C until further processing. We selected samples from collection dates that had at least 10 individual wingless adult R.padi for further testing. Similar to the samples from Kansas, DNA and RNA were extracted from individual aphids (10 individuals/field site/year) using the AllPrep DNA/RNA Mini Kit and stored at −80°C until further processing.

We generated two unique data sets for developing our applied ecological modeling framework. Both data sets use field caught aphid vectors (S. avenae and R. padi) and combine virus screening with commonly used approaches for characterizing relative abundances of microbial taxa/symbionts from insect hosts. The two microbiome data types used are DNA sequenced based but estimate symbiont abundances differently, either focusing on 1) individual taxa (qPCR data set) or 2) selected taxa identified from within a broader microbial community (metabarcoding data set).

Using a hierarchical Bayesian framework we developed a model linking virus acquisition, aphid vector and microbial symbiont abundances, and average spring temperatures. Figure 2 shows results from each focal component of our model (i.e., symbiont, virus, vector), but also refer to Figure 1 for graphical overview of our approach. Specifically, our response variables include vector abundance, virus acquisition and Buchnera relative abundance and predictor variables include spring temperature and Buchnera relative abundance. A unique feature of our models is that each microbial symbiont can be either included as a response variable (Figure 2A) or predictor variable (Figures 2C, E). Each horizontal set of panels thus demonstrates the modularity of this approach and usefulness for teasing apart complex multi-tropic relationships and the role of environmental covariates. Main outputs from these models are the coefficient estimates (e.g., slopes), which provide information on magnitude and direction of effect sizes for relationships between response and predictor variables, and 95% credible intervals (CIs) that indicate level of certainty of these effects. For relationships depicted in Figures 2A, D, E we used informative priors (i.e., assumed a positive slope) based on well established relationships between temperature, Buchnera abundance and insect development.

First, Buchnera relative abundance, measured using qPCR, shows a positive association with increased average spring temperatures (Figure 2A) across our Kansas field locations. There is an overall medium positive effect size, but the large credible intervals (CIs) suggest a moderate level of uncertainty. This positive relationship is not unexpected given the assumptions of the model and small range of average temperatures during the sampling period. Next, we evaluated relationships between virus acquisition (i.e., probability of a vector being viruliferous, tested positive for BYDV), spring temperatures (Figure 2B) and Buchnera relative abundances (Figure 2C). Here we find large effects of both predictor variables with high levels of certainty based on CIs. Interestingly, mean spring temperature had a positive coefficient estimate, while Buchnera relative abundance had a coefficient estimate that was negative. These results predict that increases in temperature during the sampling period increase the likelihood S. avenae is viruliferous, while vectors with higher Buchnera titers may be less likely to have acquired the virus. Finally, results indicate that vector abundance during the sampling spring period increases with temperature (Figure 2D) and Buchnera relative abundance (Figure 2E). Coefficient estimates are both positive, but in the case of Buchnera relative abundance there is less certainty given the large CIs. However, there was a medium to large effect of temperature, which suggests that even slight increases in early spring temperatures may increase the number of S.avenae (i.e., ~40–50% higher vector abundance with 3°C increase).

We adapted our Bayesian modeling framework (see Sections 2.4.1 and 2.4.5) for a second data set where symbiont relative abundances were estimated using standard community-level metabarcoding techniques (i.e., raw sequence counts per ASV). We expanded our approach to include two bacterial symbionts (Buchnera and Serratia) using a second aphid vector species (R. padi). Relative abundances of these symbionts varied considerably, with Buchnera being ~1000-fold higher than Serratia based on read counts. This data set differs from the one collected in KS in that there is no misalignment of data (i.e., measures of symbiont abundance and virus incidence are collocated), thus we did not utilize data fusion techniques. However, in this case we do not have information on vector abundances. Figure 3 shows results from each focal component of our model (i.e., symbiont, virus), which links virus acquisition, microbial symbiont abundances, and average spring temperatures.

Interestingly, the two vector-associated symbionts show opposite relationships with temperature during the sampling period (Figures 3A, B). Buchnera relative abundance shows a strong negative association with increased average spring temperatures (Figure 3A) across our IN field locations. In contrast, Serratia abundances are strongly positively associated with increasing temperature (Figure 3B). In both cases these are large effects sizes with moderate to high levels of certainty based on CIs. However, neither symbiont (Figures 3C, D) nor temperature (Figure 3E) was found to drive differences in the probability that R. padi was viruliferous.

How can researchers apply our modeling approach to understand microbiome mediated spread of insect-borne plant pathogens? One strength of our approach is the ability to use existing data sets that may be spatially and/or temporally disconnected. Specifically, data integration or fusion techniques allow data sets collected at different times, locations and sources (i.e., plant, insect and microbe metrics) to be combined into a single model and thus avoids the need for additional, potentially cost-prohibitive, data collection (Schaub and Kéry, 2021; Cook et al., 2023; Frost et al., 2023). Researchers can therefore take advantage of previously collected data sets that are not spatially or temporally aligned, or may be partially incomplete, to generate a system-level understanding by sharing information across data sets. In addition to using existing data, our modeling approach is informative for experimental design and directing additional data collection. This opens exiting new avenues for incorporating multiple types of microbiome data, such as gene copy numbers or absolute abundances, rapid sensor based identification (i.e., presence/absence) and even microbial metabolite or other omics data sources. The modeling tools we develop here also easily facilitate interdisciplinary research approaches by melding organisms across kingdoms with a wealth of environmental variables of interest within agroecosystems. For example, models could also begin to tease apart interactions between different insects (e.g., presence of natural enemies or competitor species) associated with the same host plant, which could in turn influence vector–microbe–plant interactions and disease dynamics.

We created a broadly applicable framework for modeling the distribution of individual microbes to identify environmental covariates (e.g., weather, landcover) that predict differences in abundance within vector microbiomes, that are then linked back to aspects of vector competence (e.g., acquisition, transmission). There are several key considerations we suggest are important for researchers when building joint vector–microbiome models using the Bayesian techniques and data integration illustrated in this study. First, identifying relevant components of vector competence specific to the focal pathosystem is necessary, particularly where there is potential for involvement of microbial symbionts (e.g., host plant preference, vector fitness). Importantly, model assumptions or inclusion of informative priors and selection of environmental covariates will depend on existing biological knowledge within each system. Bayesian modeling approaches are also able to include multiple environmental predictor variables (i.e., major and minor drivers), but researchers should be aware of potential issues with collinearity of covariates. Second, researchers should consider the level of detail or type of information provided by available data sets and when designing experiments. Our modeling approach can accommodate a range of data types that provide varying levels of detail or information about the system. For example, finer scale measurement of viral titer levels or amplicon sequencing can provide additional information beyond what is possible with dichotomous (i.e., presence/absence) viral data, such as patterns in abundances of specific viral sequence variants and titer levels within vectors. However, when addressing data deficient scenarios there are several options available, including model based solutions such as data fusion and use of informative priors based on expert knowledge or the potential need for additional data collection. For example, variation in disease incidence across growing regions and abundance of symbionts are important factors to consider for model building. In this study we demonstrate the usefulness of Bayesian techniques for modeling data sets that are relatively information rich (high disease – Great Plains region) or information deficient (low disease – Midwest region; some very low abundance symbionts, no vector data).

Here we use the aphid-BYD system as an example for developing joint vector–microbiome modeling approaches that can be customized to different pathosystems of interest. For example, a wealth of information exists on the microbial symbionts and plant pathogens vectored by whiteflies that could be used for developing Bayesian hierarchical models to investigate linkages between symbionts, environmental variables and disease dynamics (Andreason et al., 2020; Wu et al., 2022). Researchers could develop ecological models that incorporate candidate symbionts (i.e., response variable = symbiont abundance) predicted to be involved in pathogen transmission (e.g., Hamiltonella, Rickettsia) or that influence vector biology more generally (e.g., obligate symbiont Portiera). Models could be developed to compare virus and symbiont abundances across whitefly biotypes or populations with varying vector competence and could include environmental predictor variables (e.g., temperature, host plant species) that may impact virus transmission and/or acquisition. Interestingly, Rickettsia has been shown to influence whitefly thermal tolerance via increased expression of stress responsive genes (Brumin et al., 2011), making this facultative symbiont an interesting target for modeling interactions with begomoviruses. Documented variation in location of facultative symbionts within whiteflies also offers opportunities to model microbe–microbe interactions, as some symbionts are co-localized with the obligate symbiont inside bacteriocytes whiles others infect multiple organs (e.g., gut, salivary glands) and can circulate in the hemolymph (Andreason et al., 2020).

Despite extensive research on Luteovirus–aphid interactions, it remains unknown to what extent symbiotic microbes impact BYDV transmission. Early work in R. padi and S. avenae found a molecular chaperone (GroEL) produced by the primary symbiont Buchnera with high binding affinity for BYDV particles in vitro (Filichkin et al., 1997). Thus, it was proposed GroEL protects BYDV from aphid immune defenses during circulative transport through the hemolymph. However, the role of GroEL in BYDV transmission has been questioned. Bouvaine et al. (2011) demonstrated GroEL protein is restricted to specialized cells called bacteriocytes and not present in aphid hemolymph, suggesting GroEL–BYDV interactions are unlikely. More recently, studies have found Buchnera titers are lower in viruliferous aphids involved in transmission of Potato leafroll virus (Patton et al., 2021) and in aphids feeding on virus-infected soybean tissues (Cassone et al., 2015). General disruption of aphid symbiont communities using antibiotics is also linked to reduced transmission efficiency of BYDV (Yu et al., 2022). However, the role of symbionts in transmission of aphid-borne pathogens has primarily focused on Buchnera, while additional facultative symbionts have largely been ignored (Angelella et al., 2018).

Our modeling framework adds layers of mechanistic understanding to BYDV pathosystem dynamics, specifically by investigating environmental mediators of vector symbiont abundances and links to virus acquisition (Figures 2, 3). Results from our previous spatio-temporal modeling work indicated the probability an S. avenae was viruliferous (+BYDV) was unaffected by precipitation but was positively correlated with fall temperatures and distance to forest or shrubland (Enders et al., 2018). Here we further add a layer of mechanistic understanding by showing that titer levels of the obligate symbiont Buchnera are strongly negatively associated with virus acquisition (Figure 2C), and that spring temperatures may also affect Buchnera within the S. avenae aphid vector (Figure 2A) as well as probability of virus acquisition (Figure 2B). Taken together these results suggest that temperature could influence S. avenae vector competence via effects on the obligate symbiont. In contrast, R. padi virus acquisition was not linked to temperature or symbiont abundances (Figures 3C–E), suggesting the role of microbes and environmental factors are likely vector species specific. However, further in-depth laboratory studies are needed to identify exact molecular or physiological mechanisms underlying potential virus–symbiont interactions or that could explain observed reductions in Buchnera titer with S. avenae virus acquisition. Finer scale measurement of viral and symbiont titers using qPCR and investigation of variation across viral strains using amplicon sequencing could also generate additional insights from ecological models.

Our results also highlight the complexity of interactions between temperature and insect symbionts. Associations between relative abundance and temperature varied across Buchnera and Serratia (Figures 2A, 3A, B), suggesting increasing temperature may be stressful for some symbionts but not for others. Alternatively, in the case of free-living symbionts like Serratia, warmer temperatures may increase overall abundances in the environment, which could lead to higher exposure and acquisition by aphid hosts. In line with our results, a recent survey across aphid species found Serratia symbiotica prevalence was significantly greater when temperatures were higher at time of sampling (Pons et al., 2022).

Natural next steps towards a systems-level approach are to incorporate measures of plant health (e.g., disease severity, yield) and management practices (e.g., insecticides, resistant crop cultivars) into our models that are relevant for cereal agroecosystems. To develop a more complete understanding of microbial-mediated vector competence, models will also need to include BYDV transmission to new host plants in addition to vector acquisition. Fine tuning model assumptions and environmental variable selection will continue to advance our understanding of vector–BYDV–symbiont interactions. The Bayesian modeling framework applied here is flexible and can easily accommodate additional environmental covariates. For example, precipitation levels and landcover estimates (e.g., proximity to natural vs. cultivated landcover) can be easily added to existing single covariate models. Precipitation is associated with increased early season activities of some aphid vectors (Crossley et al., 2022) and is predicted to indirectly benefit vector populations by generally supporting host plant health, including growth of both cereal crops and additional grass hosts that provide a “green bridge” for vector survival between growing seasons (Thackray et al., 2009; Rashidi et al., 2021). Relatedly, exposure to a greater diversity of host plants in a given landscape could in turn contribute to increases in symbiont diversity within vectors populations (Hansen and Moran, 2014; Blankenchip et al., 2018; Enders et al., 2022). Presence of additional grass hosts are also problematic because they can contribute to viral reservoirs in the landscape that then facilitate transfer to cereal crops (Rashidi et al., 2021). Landcover (e.g., % cropland vs. non-cropland) is therefore another important candidate environmental covariate linked to the potential for aphids to acquire symbionts from alternative grass hosts that could influence vector competency and the spread of BYDV.