Data collection
Study sites
We selected four study sites in the Upper Gunnison Basin in Colorado, USA. As we intended for these sites to capture a large amount of variation in environmental conditions, the primary criteria we used for site selection was elevation. Mirroring patterns of mountain weather observed worldwide, elevation is correlated with a multitude of environmental factors in the Rocky Mountains. Predominant trends include a negative correlation between altitude and temperature, a positive correlation between altitude and precipitation34, and transitions in plant community composition (from ‘mountain shrub’ to ‘montane’ to ‘subalpine’ to ‘alpine’) with increasing altitude35.
Our sites span a large elevation gradient of approximately 1000 m that captures significant variation in environmental conditions. These sites should not be interpreted as a direct ‘elevation for climate substitution’, especially considering the significant changes in weather at each site that occur over the course of our observation periods, which ran from June to July. The lowest site, ‘cement creek’ (CC; approximately 38.82156° N, 106.86893° W), falls within a sage brush meadow on the boundary between the mountain shrub and montane vegetation zones at 2440 m elevation. The second lowest site, ‘bus turnaround’ (BT; approximately 38.97130° N, 106.99595° W), is situated in an open subalpine meadow at the base of Gothic Mountain at ~ 2940 m. The second highest site, ‘gothic mountain’ (GM; approximately 38.97969° N, 107.01937° W), is also in the subalpine zone on a steep hillside within a clearing in an evergreen forest on the lower slopes of Gothic Mountain at 3220 m. The highest site, ‘high meadow’ (HM; approximately 38.96779° N, 107.02184° W), is on an exposed meadow at the upper fringes of the subalpine zone on the shoulder of Gothic Mountain at 3,410 m.
We began our observations at each site by mapping the flax population (see “Population mapping” sub-section below). This occurred on 6/15/2020 for CC, 6/17/2020 for BT, 6/16/2020 for GM, and 6/18/2020 for HM. We concluded all observations at CC on 7/27/2020, at BT on 7/29/2020, at GM on 7/28/2020, and at HM on 7/10/2020 (aside from spore trap collection which concluded on 7/21/2020).
Weather monitoring
We deployed environmental sensors to record longitudinal weather data at each of our sites. We recorded temperature, humidity, rainfall, and wind speed and direction (at ~ 1 m above ground) every five minutes (see Supplementary Text 1 for logger specifications). We used temperature and humidity data to calculate absolute humidity (Supplementary Text 2). Due to logistical constraints, rainfall and wind speed/direction loggers were deployed later than temperature and humidity loggers.
We extracted weather metrics for a given observation period from this data as follows: We calculated the mean, maximum, and minimum temperature as the mean, maximum, or minimum of temperature readings. Likewise, we calculated mean absolute humidity as the mean of absolute humidity records (calculated from temperature and relative humidity every five minutes). Daily mean rainfall was calculated as of mean rainfall records multiplied by the number of 5 min increments in a day (i.e. 288).
Future climate data
To enable forecasting of the effects of climate change on within and between host transmission processes, we compiled projected future weather data. We extracted downscaled climate and hydrological data generated using the Localized Constructed Analogs (LOCA) method and the CESM1(CAM5) model from the “Downscaled CMIP3 and CMIP5 Climate and Hydrology Projections” archive36,37,38,39,40. Specifically, we extracted values for daily mean temperature, relative humidity, and rainfall from the CESM1(CAM5) hydrology model, and values for daily maximum and minimum temperature from the CESM1(CAM5) climate model for both the RCP4.5 and RCP8.5 emissions scenarios41. Temporally, this data spanned June through July in 2020, 2030, 2040, 2045, 2050, 2060, and 2070. Spatially, this data represents projections for the area spanning (38.9375, 39.0) latitude and (− 107.0, − 106.9375) longitude. This area is the finest grained spatial extent including the BT site for which projections were available.
Population mapping
In order to uncover drivers of population level epidemiological dynamics, it is necessary to characterize the ‘landscape of susceptibility’ in the host population. For spatially structured infectious disease systems such as flax rust, this involves documenting not only the infection status of each individual and any likely covariates of contact rate or susceptibility (e.g. size), but also the spatial arrangement of individuals. To accomplish this, we established 10 by 20 m transects at each of our field sites. We recorded the coordinates of each transect corner and the compass bearings associated with the transect axes (so that wind direction could be translated to the transect coordinate system) using a Garmin GPS unit (part #: 010-01735-10). At the beginning of the field season, we mapped the location, height, and infection status (healthy or diseased) of each plant within each transect. Plants with heights less than five centimeters were denoted as seedlings. While other metrics of size such as stem count could be used to form a more complete picture of plant size, we chose to use height alone as a proxy for size as it was the only measure that we could feasibly record for the many hundreds of plants within each transect. At the CC site, we were unable to map all healthy plants in certain areas of the transect due to time constraints, but we did record data for all diseased plants in these regions.
Epidemiological surveys
After documenting the initial epidemiological conditions, we tracked the spatiotemporal spread of flax rust within each transect. Approximately once per week, we visually inspected each flax individual in the transect for signs of disease, recording the location of any newly infected plants so they could be matched to a previously uninfected plant. We also recorded the height of all newly infected plants. We conducted these epidemiological surveys at least seven times at each site. Plants identified as infected in the initial population mapping and in these subsequent epidemiological surveys were marked with flags to ensure that they would not be repeatedly recorded as newly infected. Although we did not map the entirety of the CC transect, we covered the entire 10 by 20 m area in all epidemiological surveys.
Diseased focal plants
Within each transect, we designated a subset of diseased plants as ‘focal diseased plants’ and used them to make detailed measurements of within-host disease spread processes. These processes included changes in infection intensity, which we measured at the plant scale, as well as pustule growth, which we measured at the sub-plant scale using changes in pustule area and as a proxy.
All initially infected plants were designated as diseased focal plants, and we continued to give this designation to newly diseased plants (except in some cases, diseased plants with height of approximately 5 cm or less, as these plants were too delicate to be manipulated) until we obtained at least 25 infected focal diseased plants per site. This occurred on 7/20/2020 at CC, on 7/1/2020 at BT, and on 7/9/2020 at GM and HM. We inserted metal tags into the ground at the base of diseased focal plants so they could be quickly and reliably identified. To facilitate longitudinal measurements of disease spread at sub-plant scales, we marked up to three infected stems on each focal diseased plant with pieces of colored flagging tape tied around the base of the stem. These stems were chosen haphazardly. On each of these stems (when possible), we marked the tip of one infected leaf with black ink. We preferentially marked leaves with one or a few pustules. Approximately once every three days, we recorded detailed measurements of plant height and within-host disease spread processes.
To measure pustule growth, we photographed each marked leaf with a millimeter ruler in frame using a Cannon EOS Rebel T7i DSLR camera fitted with a Canon EF-S 35 mm f/2.8 Macro IS STM lens and a Cannon MR-14EX II Macro Ring Lite. While achieving consistent photographic angles and lighting conditions is impossible in a field setting, we attempted to keep the leaf and the scale in the same plane while maintaining a perpendicular shooting angle. We also adjusted camera settings so that pustule boundaries could be clearly distinguished in images. We used ImageJ software42 to measure the maximum and minimum diameters of each photographed pustule, and calculated pustule area assuming pustules to be elliptical. To enable the size of individual pustules to be tracked across observations, each measured pustule in each image was labeled so that it could be re-identified in subsequent images based on its location on the leaf and position relative to other pustules. We omitted data for any pustule that could not be confidently identified or whose borders became indistinct from other pustules. We marked a new leaf when we were unable to find the previously marked leaf, when the leaf was accidentally removed during data collection, or when the condition of the leaf deteriorated significantly. We suspect that either the presence of ink or our handling of the plant contributed significantly to the rate at which leaf condition deteriorated, but nevertheless we were able to observe the growth of many pustules over many weeks.
To measure infection intensity, we first counted the number of infected stems. Next, for each marked stem, we recorded the length of infected tissue as the length of the stem segment spanning the lowest (relative to the ground) pustule to the highest pustule. Finally, we counted the number of pustules present on a haphazardly selected leaf in the middle of the region of infected tissue on each marked stem. We did not use the same leaf at each observation as the middle of the region of infected tissue changed due to disease spread. Using these measurements, we calculated an infection intensity metric as the product of the number of infected stems, the average length of infected tissue, and the average number of pustules present on a leaf in the middle of the region of infected tissue.
Healthy focal plants
As we wished to investigate patterns of growth in both healthy and infected plants, we designated certain uninfected plants as healthy focal plants, and measured their height whenever we measured the height of focal diseased plants. To control for potential effects of stem marking on the growth of diseased focal plants, we made these same markings on healthy focal plants. Healthy focal plants were selected haphazardly, and we added new healthy focal plants over time to keep the number of diseased and healthy focal plants approximately equal. At some sites, the number of healthy focal plants fell towards the end of the observation period as many healthy focal plants became infected. When this occurred, we stopped recording data for the plant, and in some cases designated it as a diseased focal plant.
Spore deposition
To measure the distribution of spore deposition from an infected plant, and how it relates to the infection intensity of the source plant and wind patterns, we deployed spore traps in arrays around a subset of focal diseased plants. We chose focal disease plants that were as removed as possible from other diseased plants to minimize the number of spores originating from other plants that would be caught in the spore traps. These spore traps consisted of a ~ 2 cm2 section of Scotch Permanent Clear Mounting Tape (part #: MT76272-5) affixed to plastic backing with double sided tape. These spore traps were secured into the ground with ~ 5 cm nails. We deployed these traps at distances of either {~ 5 cm} (traps at ~ 5 cm were placed as close to the plant as possible), {~ 5 and 25 cm}, {~ 5, 25, and 50 cm}, or {~ 5, 25, 50, and 100 cm} in each of the four directions corresponding to the axes of the transect in which they were located. The traps were left in the field for approximately one week, and then collected. The sections of mounting tape were transferred to microscope slides and sealed in place using clear packing tape. We then used a light microscope to count the number of spores present in a known area. We counted spores in an area of 0.8 cm2 for most spore traps that contained few or no spores, but we counted a minimum area of 0.1 cm2 for all spore traps, including those with greater than 10,000 spores per cm2. Melampsora lini spores were identified visually using a slide prepared with spores as a reference for spore size, shape, and color. The spores matched descriptions in the literature43. We did not observe any rust fungi other than M. lini in or around our transects, so we can be reasonably certain that our counts represent only M. lini spores from infected focal diseased plants.
Statistical analyses
Plant growth and within-host disease spread
To infer how different weather factors affect plant growth, pustule growth, and changes in infection intensity, we fit generalized additive models18. We used longitudinal observations of plant height in healthy and diseased focal plants to fit the model of plant growth, and longitudinal observations of within-host disease spread at various scales in focal diseased plants to fit all other models. Because the observation periods associated with our measurements of these processes varied, we formatted the response variables in change-per-day units: For our analyses of plant growth, pustule growth, and infection intensity progression we used change in plant height per day, change in pustule area per day, and change in infection intensity per day respectively as the response variables.
To infer the effects of weather variables, we included smooth terms for these factors as predictors in each model. For each observation of change in plant height, pustule area, or infection intensity, we determined the start and end timepoints of the observation period to extract the corresponding weather data. For observations of plant height, we used 12:00:00 on the day of the observation as the start/end timepoints of the observation period. For observations of pustule area, we used photograph timestamps as the start and end timepoints. For observations of infection intensity, we used the average of the timestamps of photographs taken of pustules on that plant’s leaves as the start and/or end timepoints if such photographs were taken, because all data for an individual diseased focal plant was generally collected at once. If such photographs were not taken, we used the mean timestamp of all photographs taken at that plant’s site on the observation date as the start and/or end timepoint. To focus our analysis on fine-grained relationships between weather and within-host disease spread, we discarded data corresponding to observation windows of eight or more days. According to this criterion, 13 of 609 data points were discarded in the analysis of plant growth, 216 of 3535 data points were discarded in the analysis of pustule growth, and 25 of 338 data points were discarded in the analysis of infection intensity progression. Using the start and end timepoints of the observation period as bounds, we extracted the following variables from the weather data: temperature (mean, maximum, and minimum), absolute humidity, and mean daily rainfall.
To account for the nested nature of our data, we included smooth terms for study site and plant identity as random effects (accomplished by setting bs = ”re” in the s() function in mgcv) in each model. In the case of the analyses of pustule growth, plant identity refers to the plant on which the pustules were observed. We also included a term in each model that accounted for the value of the previous observation to capture any allometric effects. In the model of change in plant height per day, this term was a full tensor smooth product of last observed height and infection intensity. In the model of change in pustule area per day, this term was a smooth of last observed pustule area. In the model of change in infection intensity per day, this term was a full tensor smooth product of the base 10 logarithm of the last observed infection intensity and last observed maximum plant height.
We fit these models using the gam() function in the mgcv R package and the restricted maximum likelihood parameter estimation method18,44. We used gaussian response distributions with identity link functions. The dimensionality of all basis functions was set to the default value determined by mgcv. The one exception to this involved the marginal basis for log10 infection intensity (a part of the tensor term) in the model of change in infection intensity per day. The dimensionality of this marginal basis was increased to 15 after model diagnostics indicated that the default basis dimension was insufficient. To implement variable selection so that insignificant predictor terms would be effectively removed from the models, we used the double penalty approach of Mara and Wood45 (via the ‘select = TRUE’ option in the mgcv gam() function).
The coefficients of the fit models can be used to infer the effects of weather factors on plant growth and within host disease spread. To translate these inferences into predictions about how climate change might affect plant growth and within-host disease spread, we simulated trajectories of plant growth, pustule area, and infection intensity under various climate conditions. For each class of simulation, we started by defining a hypothetical starting state. For simulations of plant growth, this was a 10 cm tall uninfected plant (the plant was assumed to remain uninfected for the entire simulation). For simulations of pustule growth, this was a pustule with area 0.1 cm2. For simulations of infection intensity, this was a 25 cm tall plant with infection intensity 1.0. Next, we simulated 100 trajectories from the model posteriors using observed weather data for each site and a step size of seven days. In these simulations, we restricted plant height, pustule area, and infection intensity to remain greater or equal to 5 cm, 0 mm2, and 0.1, respectively. We included the random effect of site (but not plant identity) when making these predictions. After confirming that these simulation results qualitatively recapitulated observed patterns, we simulated sets of 100 trajectories for future climate data sets. We again excluded the random effect of plant identity when making these projections. We used the same set of 100 random number seeds for each set of 100 trajectory simulations.
While performing these simulations, we extracted weather variables as follows: When extracting weather variables from observed data sets, we defined the observation period as 00:00:00 on the first day in a step to 00:00:00 on the last day in a step (e.g. midnight on July 1st to midnight on July 8th). We extracted weather metrics as described in the “Weather data” subsection of the “Data collection” section above. Applying these same methods to the projected weather data is not possible, as these data sets only give only daily mean values. Instead, we calculated mean temperature as the mean of daily mean temperature values, maximum temperature as the maximum of daily maximum temperature values, and minimum temperature as the minimum of daily minimum temperature values. Similarly, we used daily mean temperature and relative humidity to calculate daily mean absolute humidity, and calculated mean absolute humidity as the mean of these values. We applied this same procedure to calculate mean rainfall.
Spore dispersal
To infer how spore deposition is related to source plant infection intensity, wind speed, and wind direction, we fit a tilted gaussian plume model22 (TGPM) to our spore deposition data. The equation specifying the TGPM is as follows:
$$begin{aligned} & dleft(I,H,s,X,Yright)=frac{I k {W}_{s}}{2 pi s {sigma }^{2} }{e}^{frac{{-Y}^{2}}{2{sigma }^{2}}-frac{{(H-frac{{W}_{s }X}{s})}^{2}}{2{sigma }^{2}}} & for quad X>0 0 & for quad Xle 0.end{aligned}$$
This equation describes the concentration of spores deposited at a given point ((d)) as a function of wind speed ((s)), the infection intensity of the source plant ((I)), a constant relating (I) to the source concentration of spores ((k)), and the coordinates of the point (X,Y), relative to the source. The coordinate system has the source at the origin, the X-axis parallel to the wind direction (with wind flowing in a positive direction), the Y-axis perpendicular to the X-axis on the plane defined by the ground (assumed to be flat), and the vertical Z-axis orthogonal to the X and Y axes. The shape of the three dimensional spore plum emanating from the source is defined by (s), along with constants specifying the falling velocity of spores (({W}_{s})), the height of the source ((H)), and the standard deviation of spore dispersion (({sigma }^{2})) in the horizontal and vertical directions. Following Levine and Okubu22, we calculate ({sigma }^{2}) as a function of a diffusion coefficient, (A), assuming that the variance in spore distribution increases linearly with time. We also assume that spores diffuse along the Y and Z axes at equal rates.
$${sigma }^{2}=frac{2AX}{s}.$$
We fit the parameters (k, {W}_{s},) and (sigma ) by minimizing the sum of squared differences between model predictions and the spore deposition data we collected from the spore traps arranged around diseased focal plants. We assumed that all spores deposited on a spore trap originated from the associated diseased focal plant, and that spore deposition occurred between noon on the day of deployment and ceased at noon two days post-deployment as the spore traps lost their adhesive properties and became saturated with dust. When predicting spore deposition using the TGPM, we set (I) equal to the infection intensity of the source diseased focal plant on the day of spore trap deployment and (H) equal to half of the maximum height of the source diseased focal plant at the day of spore trap deployment. As wind speed and direction were recorded every 5 min, we calculated total spore deposition as the sum of predicted deposition across five minute time windows. In each individual estimate, we set (s) equal to wind speed and calculated (X) and (Y) from wind direction and the location of the spore trap relative to the focal diseased plant (see Supplementary Text 3).
Transmission
Using the fitted TGPM we next sought to connect patterns of spore dispersal to transmission. The first step in this process was to align data on the spatiotemporal spread of disease with the locations of all healthy and infected plants (including non-focal plants) within each transect so that the relative locations of these could be determined (see Supplementary Text 4). Next, we built a dataset (‘corrected plant height data’) describing the height of all plants on the dates of epidemiological surveys. We preferentially used actual measurements from the initial population survey, epidemiological surveys, or focal plant measurements where they existed. For the majority of plants, height was not recorded on the date of a given epidemiological survey. In these cases, we fore- or hind-casted height from a previous or future record of plant height (whichever was closer) using the fitted generalized additive model of plant growth and observed weather data spanning the period between the date of interest and the closest record. We included the random effect of study site when making predictions, but not the random effect of plant identity as we wished to forecast heights of both focal and non-focal plants and only focal plant data was used to fit the model of plant growth.
After completing these steps, we next built a ‘transmission data set’ connecting weather conditions, predicted spore deposition, plant height, and infection occurrence across the time windows between each epidemiological survey. We included separate entries for each healthy plant (focal and non-focal) during each time window. We limited the scope of this data to the periods for each transect in which all infected plants (that were not seedlings) were tracked as diseased focal plants. Using noon on the date of one epidemiological survey as the beginning of the time window, and noon on the date of the next epidemiological survey as the end of the time window, we extracted the same weather variables as described in “Plant growth and within-host disease spread”. Plant height values were extracted from the ‘corrected plant height’ data set. We determined that a plant became infected if a plant with identical coordinates (after data alignment) was recorded as newly infected in the epidemiological survey that took place at the end of the observation period. Total spore deposition for each healthy plant was calculated as the sum of predicted spore deposition originating from each of the plants in the same transect that was infected at the beginning of the time window. To predict the spore deposition experienced by a healthy plant that originated from an individual diseased plant, we used the TGMP(.) We extracted the height of the diseased plant at the beginning of the time window and set (H) to half of this value. Likewise, we extracted the infection intensity of the diseased plant at the beginning of the time window, and set I equal to this value. For diseased seedlings that were not tracked as diseased focal plants, we set I = 0.1, as most were observed to be lightly infected. There were 14 instances of missing infection intensity data for a diseased focal plant on a certain date. In these cases, we fore- or hind-casted infection intensity for that plant using the generalized additive model of infection intensity progression fit in “Plant growth and within-host disease spread” (while including the random effect for study site, but not for plant identity), the closest observation of infection intensity for that plant, and the weather data spanning the time period from the missing observation to the closest measurement. We set the values of (k, {W}_{s},) and (sigma ) to those fit in the “spore deposition” section above. As wind speed and direction were recorded every five minutes, we summed spore deposition over 5 min windows that matched the resolution of the wind speed and direction data. For each of these windows, we predicted spore deposition after setting the values of X and Y based on the relative location of the healthy and diseased plants and the direction of the wind (see Supplementary Information).
Using this data set, we modeled infection outcomes in a generalized additive model framework. We assumed a binomial response distribution, coding the outcome of an infection as 1 and the outcome of no infection as 0. To account for variation in the duration of observation period, we used a complementary log–log link function, and included the log time as an offset in the model46. We included smooth terms for all weather variables as predictors, along with smooth terms study site and plant identity as random effects (accomplished by setting bs = ”re” in the s() function in mgcv). To infer the connection between spore deposition and infection, we also included a full tensor smooth product between the base 10 logarithm of mean total spore deposition per day and maximum plant height. When fitting the tensor, we fixed the smoothing parameter for the marginal smooth for log10 total spore deposition to (1times {10}^{10}). This effectively sets the shape of this marginal smooth to be linear, preventing any overfitting involving non-monotonic effects of spore deposition. Height was included as a factor in the tensor for two reasons. First, we calculated spore deposition for a point. In reality, the spore deposition experienced by a plant depends on the size of the ‘target’ that it presents. As such, larger plants are presumably challenged by a greater number of spores than a smaller plant for the same level of predicted spore deposition. Secondly, plant size could be correlated with quantitative resistance, because leaf age can influence pustule development20. To achieve computational feasibility, we fit this model using the bam() function in the mgcv R package and the fast restricted maximum likelihood parameter estimation method (“fREML”). We again implemented variable selection so that insignificant predictor terms would be effectively removed from the models via the double penalty approach.
To infer the effects of climate change and to illustrate the relationships between spore deposition, plant height, and infection risk, we predicted odds of infection using the estimated parameters of the fitted GAM (Fig. 5; Supplementary Fig. S14) and various weather data sets. In all of these projections, we included the random effect of site but excluded the random effect of plant identity.
Epidemiological model
Model description
To predict how the population level epidemiological dynamics of flax-rust might be affected in various climate change scenarios, we constructed a spatio-temporal epidemic model. This model tracks the infection status, infection intensity, and height of plants over time as infection spreads within the modeled population. There are three components required to simulate the model: (1) initial conditions describing the locations and starting states of all plants, (2) weather data, and (3) models of spore dispersal, transmission, infection intensity progression, and plant growth.
The initial conditions for individual simulations of the epidemic model were in all cases constructed to represent either the BT or GM study populations on the date of an early epidemiological survey. For BT, this date was chosen as the first which had complete weather data (spanning 00:000:00 to 24:00:00), which was 6/24/2020. For GM, this date was set to 7/2/2020. Weather data was available for earlier dates, but a sharp increase in prevalence occurred during the prior week, perhaps due to prior missed observations of new infections. The model consistently underestimated observed prevalence when initialized prior to this increase in prevalence. We did not use the HM study population to initialize the model as no epidemiological surveys occurred on a date with complete weather data. Likewise, we do not use the CC study population to initialize the model because not all regions of the transect were mapped (see “Population mapping” sub-section in “Data collection” above).
The construction of initial conditions to represent a given site proceeded as follows. The locations of all plants were taken from the ‘corrected plant height’ data set. Plants observed to be infected on the date of the first epidemiological survey were initialized as infected, and all other plants were initialized as uninfected. Following the first epidemiological survey, all infected plants (except for several that were approximately 5 cm or less in height) were designated diseased focal plants and their infection intensity was recorded. We used these measurements as the initial infection intensities of infected plants. We fore- or hind-casted infection intensities for any missing records using the method described in the transmission sub-section of the statistical analyses section above. Plant heights were taken from the ‘corrected plant height’ data set.
The weather data component of the model is used to predict spore dispersal, transmission, infection intensity progression, and plant growth in conjunction with the fitted statistical models of these processes. As such, it needs to describe mean weather conditions over the period spanning a simulation step. Additionally, because the TGPM is calibrated to make predictions over 5 min intervals, the weather data must also describe average wind speed and direction over five minute windows.
We used the fitted statistical model relating weather to the processes spore dispersal (TGPM), transmission, infection intensity progression, and plant growth to simulate these processes in the epidemiological model. These models contain random effects for both site and plant identity. The random effect of plant identity was excluded when making predictions about infection intensity progression and plant growth as random effects were only fit for focal plants. We included random effects of plant identity when making predictions using the transmission model because random effects were fit for all plants in the modeled populations.
Using these components, the model simulates changes in infection status, infection intensity and plant height in time steps of 7 days. The simulation for each time step proceeds as follows:
- (1)
First, the model predicts the number of spores deposited at the location of each healthy plant over the course of the entire time step. As all infected plants are assumed to be sources of spore deposition, the model predicts total spore deposition as the sum of spore deposition originating from each infected plant. The spore deposition originating from a single infected plant is calculated as the sum of predictions over 5 min windows spanning the time step. Individual predictions were generated using the TGPM with I is set to the infection intensity of the source plant at the beginning of the time step, H set to 0.5 times the height of the source plant at the beginning of the time step, s set to average wind speed, and X and Y set according to locations of the source and target plants and wind direction (see Supplementary Information). The parameters (k, {W}_{s},) and (sigma ) were set to the values fit to spore trap data.
- (2)
Next, for each healthy plant, the model infers the probability of that plant becoming infected from predicted spore deposition at that plant’s location and mean weather metrics using the transmission model. Following this, a random number in [0,1] is drawn, and if it is equal or lower to the probability of infection that plant becomes infected. All newly infected plants are set to have an infection intensity of 1.0.
- (3)
After simulating the transmission process, the model next simulates changes in infection intensity. For each plant, we simulated a prediction of the change in infection intensity per day by drawing randomly from the posterior of the infection intensity progression model. We used the infection intensity and height of that plant at the beginning of the time step along with mean weather metrics spanning the time step as the model inputs. We multiplied the predicted change in infection intensity per day by the length of the time step and added it to the plant’s infection intensity at the beginning of the time step to calculate the new infection intensity of the plant.
- (4)
Finally, we used a similar process to simulate change in plant height. For each plant, we simulated a prediction of the change in height per day by drawing randomly from the posterior of the plant growth model. Again, we used the height of that plant at the beginning of the time step along with mean weather metrics spanning the time step as model inputs. To calculate the new height of the plant, we multiplied the predicted change in height per day by the length of the time step and added this number to the plant’s height at the beginning of the time step.
There are two sources of stochasticity within the epidemiological model. The first of these sources pertains to the probability of infection. Although the predicted probability of infection for a certain combination of predictor values is always the same, we compare this probability to a random number to determine if infection occurs. The second source of stochasticity pertains to predicted changes in infection intensity and plant height. We predict rates of infection intensity change and plant growth by drawing randomly from the posterior distributions of the corresponding fitted models. These sources of stochasticity allow us to understand the variation in epidemiological model predictions via repeated simulation.
Simulation
To validate the epidemiological model, we simulated 10 epidemiological trajectories for the BT and GM sites and compared the results to our observations of flax-rust spread to ensure a qualitative match. For each of these sites, we initialized the model using data corresponding to the site as described above and used observed weather (from the same site) data to run the simulation. When extracting weather variables from observed weather data sets, we defined the observation period as 00:00:00 on the first day in a step to 00:00:00 on the last day in a step (e.g. noon on July 1st to noon on July 8th). We included the random effects of site and plant identity when making predictions from the transmission model, and the random effects of site when making predictions from the infection intensity and plant growth models.
To infer the effects of climate change of flax-rust epidemiology, we again initialized the model using data corresponding to the for BT and GM sites, and then simulated 10 epidemiological trajectories for each set of weather data (RCP4.5 and 8.5 emissions scenarios and the years 2020, 2030, 2040, 2045, 2050, 2060 and 2070). We extracted weather metrics as described in the “Weather data” subsection of the “Data collection” section above. We used observed records of wind speed and direction to simulate spore deposition using the TGPM as the projected data sets lack these measures. As before we included the random effects of site and plant identity when making predictions from the transmission model, and the random effects of site when making predictions from the infection intensity and plant growth models. For each set of 10 model simulations corresponding to a unique site and weather data set combination, we used the same 10 random number seeds.
Source: Ecology - nature.com
