Assessing the durability and efficiency of landscape-based strategies to deploy plant resistance to pathogens

Genetically-controlled plant resistance can reduce the damage caused by pathogens. However, pathogens have the ability to evolve and overcome such resistance. This often occurs quickly after resistance is deployed, resulting in significant crop losses and a continuing need to develop new resistant cultivars. To tackle this issue, several strategies have been proposed to constrain the evolution of pathogen populations and thus increase genetic resistance durability. These strategies mainly rely on varying different combinations of resistance sources across time (crop rotations) and space. The spatial scale of deployment can vary from multiple resistance sources occurring in a single cultivar (pyramiding), in different cultivars within the same field (cultivar mixtures) or in different fields (mosaics). However, experimental comparison of the efficiency (i.e. ability to reduce disease impact) and durability (i.e. ability to limit pathogen evolution and delay resistance breakdown) of landscape-scale deployment strategies presents major logistical challenges. Therefore, we developed a spatially explicit stochastic model able to assess the epidemiological and evolutionary outcomes of the four major deployment options described above, including both qualitative resistance (i.e. major genes) and quantitative resistance traits against several components of pathogen aggressiveness: infection rate, latent period duration, propagule production rate, and infectious period duration. This model, implemented in the R package landsepi, provides a new and useful tool to assess the performance of a wide range of deployment options, and helps investigate the effect of landscape, epidemiological and evolutionary parameters. This article describes the model and its parameterisation for rust diseases of cereal crops, caused by fungi of the genus Puccinia. To illustrate the model, we use it to assess the epidemiological and evolutionary potential of the combination of a major gene and different traits of quantitative resistance. The comparison of the four major deployment strategies described above will be the objective of future studies.


S1
Available data on the duration of latent and sporulation periods for rust diseases caused by fungi of the genus Puccinia. In all these studies, both latent and sporulation periods were determined on fully susceptible hosts under controlled conditions. Different methods to estimate duration of the latent period were used by various authors: the time from inoculation until the appearance of the first sporulating lesion (red cells), or until 50% of the lesions began sporulating (green cells). When the method is unspecified, the background of the cell is left blank.  Figure. Distribution of the latent period duration (A) and length of the sporulation period (B) of rust diseases caused by fungi of the genus Puccinia. Raw data were obtained from studies estimating these parameters on fully susceptible cereal cultivars (see S1 Table). Red and blue curves are Gamma distributions estimated from the data through maximum of likelihood (function fitdistr of the R package MASS, v7.3-47); dots indicate the mean latent period (red, 10.41 days) and sporulation period (blue, 24.37 days). The variance is 8.80 days for the latent period and 105.12 days for the sporulation period.

A B Pathogen dispersal
Based on the results of previous studies designed to estimate the dispersal kernels of P. lagenophorae [109] and P. striiformis [81,110], the power-law function has a good ability to predict the dispersal of rust spores, assuming it is isotropic (i.e. uniform in all directions). We consequently use this function in our model: where a > 0 is a scale parameter and b > 2 determines the weight of the dispersal tail. The expected dispersal distance is given by: . However, the estimation of parameters a and b is not straightforward. Indeed, many studies dealing with the dispersal of rust spores did not use a power-law function [78,111], or they did but used a unidimensional or modified power-law functions [81,109,110,112,113], which makes calibration of our own function difficult. We therefore focused on the results obtained in four studies [78,109,110,113] which used exponential functions as dispersal kernels. In spite of a lower goodness of fit of the exponential function to model spore dispersal [81,110,113], estimation of a mean dispersal distance from this function is easier and enables direct comparison between studies. The overall mean dispersal distance in these four studies is 19.46 m (S2 Table), thus we calibrated a and b of our powerlaw function in such a way that μexp=20 m, which represents 1% of landscape length. Among the possible solutions we arbitrarily chose a=40 and b=7 (S8 Figure). Figure. Two-dimensional representation of the power-law dispersal kernel used in this study (μexp=20 m; a=40; b=7.). The top panel indicates the logarithm of the probability to disperse from the origin to any point of the landscape, and the bottom panel indicates the cumulative probability to disperse over a given distance. Table. Mean dispersal distance of rust spores computed from available results obtained in studies which used exponential functions. Let the dispersal function be written as ( ) = . − with x the distance of dispersal. Then the mean dispersal distance is simply 0 = 2. . The parameters of the sigmoid curve for contamination of hosts by propagules (κ=5.33, σ=3) have been parameterised as in previous works [40,42,55], such as π(0)=0, π(1)=1, and the inflexion point is located at 0 = (( − 1)⁄ ) 1 ⁄ ≈ 0.5. Thus, the contamination of a healthy host is easier when the proportion of healthy hosts is higher than 50%, and harder otherwise (S9 Figure).

A B Evolutionary parameters
Very few empirical data are available from rust pathogens to help parameterise the model with respect to mutation probabilities, cost of infectivity, cost of aggressiveness and number of mutations required to completely erode a trait for quantitative resistance. It is important to remember that the mutation probability in the model refers to the probability that a spore has a different phenotype from its parental lesion. For instance, τ1 gives the probability for a spore to have an infective phenotype on a resistant cultivar carrying major gene 1, given that this spore was produced by a lesion triggered by a non-infective pathogen. This probability depends on the number of mutations per generation per base pair (i.e. the classic 'genetic mutation rate' of empirical studies), the number and nature of the specific genetic mutations required to overcome major gene 1, and the potential dependency between these mutations. In addition to the lack of empirical data, all these parameters are likely to vary depending on the pathosystem as well as the resistance sources.
Consequently, this work focuses on a simple theoretical case, where evolutionary parameters have been arbitrarily fixed. Since one of the objectives of this application case is to identify promising combinations of qualitative and quantitative resistance, and based on preliminary simulations, the mutation probabilities were set at τg=τw=10 -4 in such way that a cultivar carrying a single major gene would be overcome in less than one year. Using these values, pathogen evolution is quick enough to potentially adapt to the resistance genes, but slow enough to allow comparison of the effectiveness of different resistance combinations with regard to controlling disease. The resistance efficiency is set at ρg=1 for qualitative resistance (i.e. a major resistance gene confers complete immunity against nonadapted pathogens), and ρw=0.5 for quantitative resistance traits (i.e. infection rate, sporulation rate, duration of the sporulation period, or number of epidemic cycles in a cropping season of a nonadapted pathogen on resistant hosts is reduced by 50%). The cost of infectivity and the cost of aggressiveness were set at θg=θw=0.5, the trade-off strength at βw=1 (linear trade-off), and the number of pathotypes with regard to the aggressiveness components at Qw=6 (i.e. every mutation step improves or degrades an aggressiveness component by 20%). This last value is a compromise between the high number of steps required to simulate a gradual pathogen evolution, and having a smaller number of different pathotypes to limit the computational time required to perform the simulations.
This parameterisation gives the following infectivity and aggressiveness matrices: Fully aggressive agw(p)=6 0.5 1 Grey lines indicate pathogen population at the beginning of the simulations.
The parameterisation of the mutation probabilities (τg=10 -4 and τw=10 -4 ) and the number of steps to completely erode a quantitative resistance (given by Qw) give the following mutation matrix for infectivity gene g:

Discussion of model parameterisation and assumptions
Parameterisation to rust pathogens.
Pathogen aggressiveness. This model has been calibrated for rust-like pathogens using available knowledge of the epidemiology of fungi in the genus Puccinia. In particular, we found many data related to pathogen aggressiveness and dispersal. This allows the reliable estimates of parameters emax, γmin, γvar, Υmax and Υvar. However, as noted above, the calculation of the effective number of spores produced daily by a lesion and effectively dispersed to a leaf where they may trigger an infection (rmax in our model) is extremely challenging. Acquisition of data to inform this parameter is crucial, since this parameter has a strong influence on epidemic spread and pathogen evolution [45]. In this study, rmax directly affects pathogen population size, which influences both the number of new infected hosts from a single infectious source, and the probability that some of these new infections are triggered by mutant spores. For example, smaller values should lead to higher times to appearance of mutants (durability measure (d1)) than what we observed in our simulations (see white bars in Fig  5A).
Pathogen dispersal. The dispersal of fungal spores has been extensively studied, and we now have considerable evidence to support the use of a power-law kernel to simulate spore dispersal [81,109,110]. Available data indicate that the mean dispersal distance is of the order of 20 m (see S2 Table). However, sensitivity analyses of models simulating plant epidemics show that both the mean dispersal distance [40,114] and the width of the tail of the dispersal kernel [42,79] have a great impact on pathogen spread. In the present work, the calibration of parameters a and b gives an appropriate mean dispersal distance, but is probably less accurate with respect to simulation of long-distance dispersal events (determined by the width of the tail of the dispersal kernel). In the future, we need studies built with exactly the same equations, to better compare and calibrate parameter values of dispersal kernels. This will be particularly important in future modelling studies aiming to explore the effect of landscape composition (e.g. proportion and spatial aggregation of a resistant cultivar) on pathogen spread and evolution (see 'Future research directions' in the discussion section).
All healthy sites of a plant do not have the same propensity to be infected [78]. This is at least partly due to plant architecture, which makes healthy sites not equally accessible to spores following dispersal. Thus, once the most accessible sites are infected, the remaining healthy sites may be harder to reach and infect. Our contamination function (see S9 Figure) accounts for this process. The high flexibility of this function makes it relevant for a wide range of pathogens, including rusts. However, we did not have any experimental data to parameterise it for cereal crops. Therefore, we used the same arbitrary values (κ, σ) as in previous works [40,55] in order to facilitate comparisons between studies. Nevertheless, since this function only mitigates the overall infection rate (emax), it is likely that epidemic spread depends more on the value of the latter than on the parameters or shape of this sigmoid function. A sensitivity analysis of a similar simulation model showed that use of a linear contamination function, rather than a sigmoid function [41] has little effect on the level of pathogen adaptation to quantitative resistance.
Seasonality and initial conditions. Parameters related to host growth (δv), host density ( ), initial conditions ( 0 and ϕ), and seasonality (λ) are specific to the agricultural (e.g. choice of cultivars, planting density) and epidemiological (e.g. initial level of contamination, presence of a wild reservoir or volunteer host plants for the pathogen) contexts. In our application case, these parameters have been parameterised so as to simulate regular host dynamic and epidemics. This helps compare different resistance deployment strategies in standardised conditions, and attribute variation in model outputs to the deployment strategies rather than to the initial conditions or the epidemiological context. On the other hand, it could be interesting to vary the growth rate of resistant cultivars (or their contribution to yield) to study the impact of cost of resistance on the performance of different deployment strategies. One could also vary the off-season survival probability of infectious hosts (λ) to investigate the potential of strategies based on green bridge suppression (e.g. fungicide treatments in the end of the cropping season; removal of potential volunteer plants during the off-season).
Pathogen evolution. As discussed previously, parameterisation of evolutionary processes in the model from empirical data is extremely difficult. In addition to the scarcity of studies designed to estimate such parameters for rusts, these parameters may not correspond with the definition of the parameters in our model (e.g. our 'mutation probability' is different from the classic 'mutation rate'; see above). Consequently, we chose to investigate combinations of qualitative and quantitative resistances in a simple context, using the same mutation probabilities (τg and τw), costs of infectivity (θg) and aggressiveness (θw), trade-off strengths (βw), and numbers of steps to completely erode a quantitative resistant trait (Qw) for every major resistance gene g and trait of quantitative resistance w.
In this model, the mutation probability gives the probability of appearance of mutants. Therefore it has a great effect on the time to appearance of such mutants [30,39] or the speed of erosion of quantitative resistance [56]. However, provided that mutation probabilities are the same for all infectivity genes and aggressiveness components (i.e. = = , ∀ , ∀ ), the specific value of these mutation probabilities should similarly affect all scenarios of our numerical experimental design. In other words, they should not change the relative durability and epidemiological efficiency of these simulated scenarios, and consequently the ranking of the tested resistance combinations. To test this hypothesis, we replicated the simulations with smaller mutation probabilities: τg= τw={10 -5 ; 10 -6 ; 10 -7 }. The results of these new simulations show that the durability of the major gene increases with decreasing mutation probabilities (S5ACE Figure), but the tested resistance combinations keep the same ranking. Similarly, the combinations of qualitative and quantitative resistances result in better epidemiological outcomes when the mutation probabilities decrease (S5BDF Figure). Nevertheless, regardless the mutation probabilities, the most promising combination is still the pyramid of two major resistance genes, followed by a major gene combined with a quantitative resistance trait against the latent period, next with a quantitative resistance trait against the infection rate and the sporulation rate, and finally with a quantitative resistance trait against the sporulation duration.
In contrast, the fitness costs associated with pathogen adaptation may differentially affect the simulated scenarios. Several studies have demonstrated that higher infectivity costs increase resistance durability [29,30] and epidemiological efficiency [33]. In our results, the durability of the pyramid of two major genes is likely due to the cost of infectivity (θg) endured by mutant pathogens (see main text and Fig 5A). However, smaller values would result in a reduced performance of this strategy in comparison to other resistance combinations. This may also be true for comparisons of pyramiding with other spatio-temporal deployment strategies such as mosaics [50]. Similarly, the time to establishment of mutant pathogens in a host population carrying quantitative resistance has been shown to increase with higher costs of aggressiveness [37] or stronger trade-off relationships [41,45,55]. Thus, the associated parameters (θw and βw), as well as the number of steps to erode a quantitative resistant trait (Qw) should considerably impact model outputs like final level and speed of erosion of quantitative resistance.
To conclude, the impact of all these evolutionary parameters, related to the choice of the resistance source, on the evolutionary and epidemiological outcomes of the simulated deployment strategies need to be explicitly evaluated within a dedicated numerical experimental design. This is the objective of future studies (see section 'Future research directions' in main text). Furthermore, given the variability in fitness costs and trade-off relationships associated with the evolution of pathogens, including rust fungi [65], a careful assessment of these parameters is required before the deployment of genetic resistance in the field, in order to get reliable predictions of their durability and efficiency.

Parameterisation to necrotrophic fungal pathogens.
Parameterisation of our model is flexible enough to encompass a wide range of pathogens. As an example, although our application case focused on biotrophic fungi, necrotrophic pathogens could be simulated, by allowing transmission from stubble after crop harvest instead of from living tissue during the cropping season. To this end, the latent period could be increased to cover the whole cropping season, and the off-season survival probability could be adjusted in such a way that the total pool of produced spores is available for the following cropping season (after which the pool of spores would no longer be available). This scenario would at least partially represent major features of interactions such as that between canola and blackleg (Leptosphaeria maculans) [115].