An epi‐evolutionary model for predicting the adaptation of spore‐producing pathogens to quantitative resistance in heterogeneous environments

Abstract We have modeled the evolutionary epidemiology of spore‐producing plant pathogens in heterogeneous environments sown with several cultivars carrying quantitative resistances. The model explicitly tracks the infection‐age structure and genetic composition of the pathogen population. Each strain is characterized by pathogenicity traits determining its infection efficiency and a time‐varying sporulation curve taking into account lesion aging. We first derived a general expression of the basic reproduction number R0 for fungal pathogens in heterogeneous environments. We show that the evolutionary attractors of the model coincide with local maxima of R0 only if the infection efficiency is the same on all host types. We then studied the contribution of three basic resistance characteristics (the pathogenicity trait targeted, resistance effectiveness, and adaptation cost), in interaction with the deployment strategy (proportion of fields sown with a resistant cultivar), to (i) pathogen diversification at equilibrium and (ii) the shaping of transient dynamics from evolutionary and epidemiological perspectives. We show that quantitative resistance affecting only the sporulation curve will always lead to a monomorphic population, whereas dimorphism (i.e., pathogen diversification) can occur if resistance alters infection efficiency, notably with high adaptation costs and proportions of the resistant cultivar. Accordingly, the choice of the quantitative resistance genes operated by plant breeders is a driver of pathogen diversification. From an evolutionary perspective, the time to emergence of the evolutionary attractor best adapted to the resistant cultivar tends to be shorter when resistance affects infection efficiency than when it affects sporulation. Conversely, from an epidemiological perspective, epidemiological control is always greater when the resistance affects infection efficiency. This highlights the difficulty of defining deployment strategies for quantitative resistance simultaneously maximizing epidemiological and evolutionary outcomes.


| INTRODUC TI ON
Resistance to parasites, that is, the capacity of a host to decrease its parasite development (Raberg et al., 2009), is a widespread defense mechanism in plants. Qualitative resistance usually confers disease immunity such that the parasite disease phenotype follows a discrete distribution ("causes disease" vs. "does not cause disease") on a resistant plant (McDonald & Linde, 2002). Quantitative resistance leads to a decrease in disease severity (Poland et al., 2009;St. Clair, 2010) such that parasite pathogenicity follows a continuous distribution (Lannou, 2012;McDonald & Linde, 2002;St. Clair, 2010).
Pathogenicity, also known as aggressiveness, can be estimated in laboratory experiments through the measurement of a number of pathogenicity traits (Lannou, 2012) expressed during the basic steps in the host-pathogen interaction. Quantitative resistance has attracted attention in plant breeding as a means of pathogen control in low-input cropping systems, particularly due to its supposed higher durability than qualitative resistance (Niks et al., 2015). However, plant pathogens also adapt to quantitative resistance (see Pilet-Nayel et al. (2017) for a review). The resulting gradual "erosion" of resistance effectiveness (McDonald & Linde, 2002) corresponds, from the pathogen side, to a gradual increase in pathogenicity.
Adaptive dynamics (Dieckmann, 2002;Geritz et al., 1997Geritz et al., , 1998 approaches assume that epidemiological and evolutionary processes unfold over different timescales. They essentially focus on longterm predictions for endemic diseases. Fewer studies (Bourget et al., 2008;Iacono et al., 2012;Rimbaud et al., 2018) have addressed the fundamental short-and long-term objectives of the sustainable management of plant diseases Zhan et al., 2015): The short-term goal is the reduction in disease incidence, whereas the longer-term objective is to reduce pathogen adaptation to resistant cultivars. Evolutionary epidemiology analysis is well suited to this purpose (Day & Proulx, 2004). Essentially inspired by quantitative genetics, it accounts for the interplay between epidemiological and evolutionary dynamics over the same timescale. As such, it can be used to monitor the dynamics of epidemics and the evolution of any set of pathogen life-history traits of interest simultaneously. It can also take into account the heterogeneity of host populations resulting, for example, from differences in the genetic, physiological, or ecological states of individuals . This is typically the case with field mixtures, in which several cultivars are cultivated in the same field, and with landscape mosaics, in which cultivars are cultivated in different fields (Rimbaud et al., 2021).
In this article, we follow this approach and study the evolutionary epidemiology of spore-producing pathogens in heterogeneous agricultural environments. Plant fungal pathogens (sensu lato, i.e., including oomycetes) are typical spore-producing pathogens responsible for nearly one third of emerging plant diseases (Anderson et al., 2004).
Spore production is usually a function of the time since infection, due to lesion aging (Kolnaar & Bosch, 2001;Sache et al., 1997;Segarra et al., 2001;van den Bosch et al., 1988). We first formulated a general model explicitly tracking the infection-age structure and genetic composition of the pathogen population. Mathematically, this model is an extension to heterogeneous plant populations of an integrodifferential model introduced by Djidjou-Demasse et al. (2017). We then investigated how the deployment of quantitative resistances altering different pathogenicity traits affected pathogen population structure at equilibrium. This question was addressed by highlighting the links between the frameworks of evolutionary epidemiology and adaptive dynamics. We characterized the evolutionary attractors of the coupled epidemiological evolutionary dynamics while emphasizing the differences between the cornerstone concepts of ℛ 0 in epidemiology (Diekmann et al., 1990;van den Driessche & Watmough, 2008) and invasion fitness in evolution (Dieckmann, 2002;Diekmann et al., 2005;Geritz et al., 1998;Metz et al., 1996;Nowak & Sigmund, 2004). Finally, we investigated the effect of deploying quantitative resistances on the transient behavior of the coupled epi-evolutionary dynamics, both at epidemiological (disease incidence) and at evolutionary (resistance durability) levels.

| AN EPI -E VOLUTIONARY MODEL FOR S P ORE-PRODUCING PATHOG ENS
Let us suppose that we cultivate two host types differing in terms of quantitative resistance level to a pathogen (e.g., a resistant and a susceptible cultivar). How can we model the joint epidemiological and evolutionary dynamics of the host-pathogen interaction? In this section, we first formulate a general model for the dynamics of spore-producing plant pathogens in structured populations, which we then apply to the study of specific scenarios.

| Host and pathogen populations
Let us consider a spore-producing plant pathogen infecting a heterogeneous host population with N c host classes. In our study, the host classes are assumed to represent different plant cultivars, but, more generally, host heterogeneity may correspond to differences in host developmental or physiological states or different habitats. At a given time t, hosts in class k can be either healthy (H k (t)) or infected.

K E Y W O R D S
adaptive dynamics, basic reproduction number, integro-differential equations, quantitative resistance, resistance durability, spore-producing pathogens Note that, in keeping with the biology of fungal pathogens, we do not track individual plants, but instead focus on leaf area densities (leaf surface area per m 2 ). The leaf surface is viewed as a set of individual patches corresponding to a restricted host surface area that can be colonized by a single pathogen individual (i.e., no co-infection is possible). Thus, only indirect competition between pathogen strains for a shared resource is considered. Spores produced by all the infectious tissues are assumed to mix uniformly in the air to form a well-mixed pool of spores that can land on any host class according to the law of mass action. Thus, the probability of contact between a spore and host k is proportional to the total healthy leaf surface area of this host.
The parasite population is structured by a continuous phenotype x and by the time since infection a, such that I k (a, t, x) represents the density of infected tissue of host k at time t, which have been infected for a duration a with strain x. The density of the airborne pool of spores with phenotype x at time t is denoted by A (t, x). Both the phenotype and time since infection affect two pathogenicity traits in each host class summarizing the basic steps of the disease infection cycle: (i) infection efficiency k (x), that is, the probability that a spore deposited on a receptive host surface produces a lesion; and (ii) the sporulation curve r k (a, x). In line with the biology of plant fungi (Van den Bosch, 1988;Kolnaar & Bosch, 2001;Sache, 1997;Segarra et al., 2001), we assumea gamma sporulation curve defined by four parameters: (i) the latent period k (x), (ii) the total number of spores p k (x) produced by a lesion throughout its entire infectious period, and (iii-iv) the rate and shape parameters of the gamma distribution (denoted k (x). and n k (x) respectively). These hypotheses lead to the following sporulation function: where Γ is the gamma function. With this formalism, the duration of the infectious period can be estimated (based on a Gaussian approximation) as 4 √ n k (x)∕ k (x).

| Epi-evolutionary dynamics
With these assumptions, we can derive the following integrodifferential equations to describe the epidemiological and evolutionary dynamics of the host and pathogen populations (see Figure 1C for a schematic representation of the model with 2 hosts): Healthy hosts are produced at rate Λ and k is the proportion of the host k at planting in the environment. Healthy hosts can become infected by airborne spores. Thtotal force of infection on a host k where ℝ N is the phenotype space of dimension N. Airborne spores produced by infected hosts become nonviable at rate . Healthy hosts die at rate (regardless of their class), and infected hosts die at rate + d k (a, x), where d k (a, x) is disease-induced mortality. Hosts infected with strain y produce airborne spores with phenotype x at rate m (y − x) r k (a, y), where m (y − x) is the probability of mutation from phenotype y to phenotype x. Thus, mutations randomly displace strains in the phenotype space at each infection cycle (i.e., generation) according to a mutation kernel m . A centered multivariate Gaussian distribution with standard deviation is typically used hereafter. However, other mutation kernels can be used provided that they satisfy some general criteria (Appendix S2).
In the simplest scenarios, the phenotype space can be onedimensional (i.e., N = 1), but we can also consider a multidimensional phenotype space (i.e., N ≥ 2). We refer the reader to Table 1 for a list of the variables and parameters of the model, and to Appendix S3 for a discussion of how to recover several simpler models in the literature from model (2).

| BA S IC REPRODUC TI ON N UMB ER AND INVA S I ON FITNE SS
Model (2) allows us to track the coupled epidemiological and evolutionary dynamics of the host-pathogen interaction. We can use a numerical integration of this model to determine the transient dynamics toward potential evolutionary attractors. For further analytical progress, we can consider two limiting cases of this process.
First, in an initially uninfected population, whether or not a single pathogen strain spreads in the population can be determined by calculating the basic reproduction number of this strain. Second, the long-term evolutionary endpoints can be determined by calculating the invasion fitness of a rare mutant in a resident population at equilibrium according to the standard adaptive dynamics methodology.

| Basic reproduction number: Invasion in an uninfected population
The basic reproduction number, usually denoted ℛ 0 , is defined as the total number of infections arising from one newly infected individual introduced into a healthy (disease-free) host population (Anderson, 1991;Diekmann et al., 1990). It is typically used to study the spread of a pathogen strain x in an uninfected host population.
In an environment with N c host classes in which all the pathogen propagules pass through a common pool in compartment A as in model (2), ℛ 0 is the sum of pathogen reproduction numbers for each host class (Rueffler & Metz, 2013). A pathogen with phenotype x will (2) m (y − x) r k (a, y) I k (t, a, y) da dy.
F I G U R E 1 (a) Shapes of gamma sporulation function r k (a, x) with latent period k = 10 and total spore production over the entire infectious period p k = 5.94 as a function of n k and k . The infectious period is ≅ 24 days for both functions. (b) Phenotypic landscapes of a pathogenicity trait on a susceptible (S) and on a resistant (R) cultivar. Trait values are described by unnormalized Gaussian functions. The optimal parasite phenotypes S and R have maximum tras max and max 1 − Δ max on the S and R cultivars, respectively. The relative effectiveness of resistance R eff = S S − R S ∕ max corresponds to the trait difference between the cultivars of phenotype S . The relative cost of adaptation C adp = S S − S R ∕ max corresponds to the trait difference between the phenotypes S and R on the S cultivar. Here, R eff = 0.8, C adp = 0.9, and Δ max = 0.2. C Schematic representation of the model with 2 hosts (a susceptible (S) and a resistant (R) cultivar) and a one-dimensional mutation kernel (i.e., N c = 2 and N = 1). All state variables and parameters are listed in Table 1 ( where the quantity ℛ k 0 (x) = Ψ k (x) H 0 k is the basic reproduction number of a pathogen with phenotype x in host k. This expression depends on the disease-free equilibrium density of healthy hosts H 0 k = k Λ∕ and Ψ k (x), the reproductive value of a pathogen with phenotype x landing on host k (i.e., the fitness function, Appendix S4). It is given by In this expression, multiplying the probability that a lesion is viable at time since infection a by r k (a, x) and integrating over all times since infection a gives the total number of spores effectively produced by a lesion during its lifetime. This differs from (x), which is the number of spores potentially produced if the lesion remains viable during the whole infectious period. A closely related expression with a simpler model (single host type and no disease-induced mortality) was derived by Gilchrist et al. (2006).
Assuming that d k does not depend on the time since infection and a gamma sporulation function (1), we obtain Furthermore, if r k is a constant (i.e., r k (a, x) = p k (x) and d k (a, x) = d k (x) for all a), we recover the classical expression of ℛ 0 .
for SIR models (Day, 2002) with

| Invasion fitness and long-term evolutionary attractors
Once a pathogen strain x has reached an endemic equilibrium, invasion by a new mutant strain can be studied with the adaptive dynamics methodology. A rare mutant strain with phenotype y will invade a resident pathogen population with phenotype x if its invasion fitness f x (y) > 0. Potential evolutionary endpoints can be identified as the solutions of f x (y) > 0. The sign of this two-dimensional function is classically visualized on a pairwise invasibility plot (PIP) (Dieckmann, 2002;Diekmann et al., 2005;Geritz et al., 1997Geritz et al., , 1998Metz et al., 1996;Nowak & Sigmund, 2004). In our model, the existence of a common pool of pathogen propagules allows us to write (Appendix S6) the invasion fitness f x (y) as Environmental feedback of the resident strain x conditions the ability of a mutant strain y to invade the resident population. It depends on the conditions established by the resident strain, particularly as concerns the resources of host k already taken by x. When infection efficiencies do not differ between host classes (i.e., k = , for all host classes k), this feedback term is unique and comes out of the sum, and using equation (3), equation (6) can be rewritten as It follows that model (2) admits an optimization principle based on ℛ 0 (Gyllenberg, & Service, 2011;Lion & Metz, 2018;Mylius & Diekmann, 1995;Metz et al., 2008;Rueffler & Metz, 2013). Indeed, the sign of the invasion fitness f x (y) is given by the sign of the difference between ℛ 0 (y) and ℛ 0 (x); and thus, the evolutionary attractors of model (2)  Accordingly, if some plant hosts affect infection efficiency, the calculation of invasion fitness with equation (6) is required for the characterization of evolutionary attractors.

| C A S E S TUDY: DEPLOYMENT OF QUANTITATIVE RE S IS TAN CE S
As an application of our general model, we now consider two habitats corresponding to a susceptible (S) and a resistant (R) cultivar differing by a single quantitative resistance trait. We assume that, after a long period of monoculture of the S cultivar, a fraction of the S cultivar is replaced by the R cultivar at t = 0. We consider two scenarios, depending on how the resistance trait affects the life cycle of the pathogen. In the SP scenario, the resistance affects only total spore production, whereas in the IE scenario, the resistance affects only infection efficiency. The analysis of the scenarios makes use of analytical results from adaptive dynamics and simulation results from evolutionary epidemiology while highlighting the links between these frameworks.

| Scenario simulations and model outputs
Phenotypic landscapes on the S and R cultivars. Model simulation first requires the specification of fitness functions describing the pathogenicity trait values of any strain x on each cultivar ( Figure 1B). For a trait (either total spore production in the SP scenario, or infection efficiency in the IE scenario), we used the unnormalized Gaussian function prob.that a lesion is viable at age a da, .
Without loss of generality, and to simplify the presentation, the optimal phenotypes are opposite for the S and R cultivars.
Therefore, the optimal parasite phenotypes S = − and R = are characterized by maximal trait values max and max 1 − Δ max on the S and R cultivars, respectively, where Δ max is the relative difference in maximal traits between the two cultivars. The inverse of the variances 1∕ 2 S and 1∕ 2 R defines the selectivity of each habitat. Rather than using these parameters, the phenotypic landscapes can be reparameterized to fit the terminology used in plant pathology to describe quantitative resistances with the relative effectiveness of resistance R eff and the relative cost of adaptation C adp ( Figure 1B).
More precisely, for a given S > 0, R eff and C adp both in 0, 1 and Δ max in 0, R eff , we have Parameter values and initial conditions. For each scenario, simulations were used to explore how the relative effectiveness of resistance R eff (between 0.5 and 0.99), the relative cost of adaptation C adp (between 0.3 and 0.99), and the deployment strategy (between 0.05 and 0.95) affect coupled epidemiological and evolutionary dynamics. Note that, within the parameter space considered, the tradeoff function that links S and R in both cultivars can be concave, convex, or sigmoidal ( Figure S1). Similarly, the basic reproduction ratio ℛ 0 can have either one global maximum, or one global and one local maximum, and this global maximum can be closer to S or to R (Figure 4, line 3). These features strongly affect the dynamics. Simulations were initiated with a density of airborne spores A (0, x) at mutation-selection equilibrium in an environment in which only the S cultivar was grown.
All model parameters and initial conditions are summarized in Table 2. Parameters for pathogenicity traits do not fit a particular pathogen species but are instead typical of biotrophic foliar fungal diseases, such as wheat rusts on susceptible cultivars. Based on Rimbaud et al. (2018), we set the infection efficiency to 0.2 and the duration of the latent and infectious periods to 10 and 24 days, respectively. Total spore production was set so as to obtain an ℛ 0 of 30 in an environment containing only the S cultivar (Mikaberidze et al., 2016).
Evolutionary and epidemiological outputs. The evolutionary ouut considered was the time to emergence T emg of the adapted strain * . This corresponds to the time from which its proportion in the airborne pool of spores relative to S remains ≥0.05 (Table 2). If the equilibrium is monomorphic, the adapted strain * is the only evutionary attractor. If the equilibrium is polymorphic, the adapted strain * is considered as the one closest to R . The epidemiological control provided by deployment of the R cultivar was assessed by determining the relative healthy area duration (rHAD) gain, a variable used as a proxy for crop yield (Iacono et al., 2012;Rimbaud et al., 2018). It was calculated over 750 generations by assuming the deployment of the R cultivar over a period of 50 years, with 15 generations of the pathogen per year (Table 2).

| Evolutionary outcomes
In the SP scenario, the only pathogenicity trait targeted by the quantitative resistance is total spore production. As the ℛ 0 optimization principle holds, the population always becomes monomorphic around an evolutionary attractor * corresponding to the global maximum of the ℛ 0 function. This happens regardless of the shape of ℛ 0 , whether ℛ 0 has a single maximum (Figure 2, line 1) or an additional local maximum (Figure 2, line 2).
TA B L E 1 Main notations, state variables, and parameters of the model. The general model is structured by time since infection a, pathogen strain x, and host class k. Some of its parameters are functions of these structuring classes

Category Description Unit
Notations In the IE scenario, the only pathogenicity trait targeted by the quantitative resistance is infection efficiency. The optimization principle based on ℛ 0 does not hold here, and evolutionary branching and diversification are possible and cannot be characterized by the shape of ℛ 0 . Typically, a polymorphic pathogen population can be selected at equilibrium, while ℛ 0 has a single maximum ( Figure S2).
In this scenario, we need to calculate the invasion fitness to characterize the evolutionary attractors (equation 6).
As in the SP scenario, the population can evolve to a monomorphic equilibrium (Figure 2, line 3). This outcome occurs for a broad range of resistance genes and deployment strategies (Figure 3).
However, a dimorphic population can also be selected (Figure 2, TA B L E 2 Initial conditions and parameters used for simulations with the S and R cultivars. Two scenarios were considered: The R cultivar can alter either the infection efficiency of the pathogen (IE) or its total spore production (SP Relative healthy area duration gain ∫ b max = max for the IE scenario and max = p max for the SP scenario. The resistance impacts only spore production. The fitness function is unimodal (left panel). At t = 0, the pathogen population is at its mutation-selection equilibrium on the S cultivar (light green distribution). The infection dynamics and phenotypic composition of the pathogen population on the S and R cultivars are displayed in the central and right panels, respectively. The black line in the right panel corresponds to the time to emergence T emg . Line 2. As in line 1, but with R eff = 0.99. The fitness function is bimodal, with both a global and a local maximum. Line 3. The resistance affects only infection efficiency. The pairwise invasibility plot (PIP) visualiz the sign of invasion fitness f x (y) (left panel). As a mutant strain y will invade the resident population x only if f x (y) > 0 , the PIP reveals a single evolutionary attractor * (the vertical line through * . is completely contained within a region marked "-"). Line 4. As in line 3 but with R eff = 0.99. The PIP reveals two evolutionary attractors (left panel). For all panels, all the other parameters are set to their reference values ( Table 2). Note that the timescale axis varies between lines number of production situations for the R cultivar characterized by the highest cost of adaptation C adp = 0.99 . In this case, the deployment of a proportion of the R cultivar above a threshold that increases with decreasing resistance effectiveness generally leads to the selection of a dimorphic population (Figure 3 . This is the case mostly when the evolutionary attractor is closer to S (optimal phenotype on the S cultivar) than to R (optimal phenotype on the R cultivar). In such settings, unless it already exists from the beginning in the mutationselection equilibrium, the evolutionary attractor can be rapidly at-

| Epidemiological outcomes
Long-term epidemiological control is assessed by determining the relative increase in healthy area duration. In the SP scenario, it ranges from 1.25 to 3.21, with a mean of 1.61 over the whole parameter space explored ( Figure 5, line 1). Achieving a substantial relative yield gain ≥1.5 requires the deployment of an R cultivar with a cost of adaptation ≥0.9 in at least 50% of the fields, regardless of resistance effectiveness. The achievement of a higher relative yield gain ≥2 requires, in addition, a high resistance effectiveness ≥0.9.
In the IE scenario, relative yield gains range from 1.27 to 5.02, with a mean of 1.9 ( Figure 5, line 1). In all the production situations explored, relative yield gains were higher for the IE scenario than for the SP scenario. The mean difference in relative yield gains between production situations was 0.28 (maximal difference 1.93).
Substantial gains ≥2 are obtained for larger ranges of resistance effectiveness and proportions of the R cultivar at planting, provided that the costs of adaptation are ≥0.9. Remarkably, for both scenarios, epidemiological control is not correlated with time to emergence in the production situations explored (compare Figures   4 and 5). In our case study, we show that the evolutionary and epidemiological F I G U R E 4 Time to emergence of the evolutionary attractor. The time to emergence T emg is the time at which the proportion of the evolutionary attractor (the closest to the optimal phenotype on the R cultivar) becomes ≥5% in the pool of spores in the air. Line 1. The resistance affects only total spore production (SP scenario). T emg is characterized as a function of the proportion of the R cultivar at planting (x-axis) and of the relative effectiveness of the resistant cultivar (y-axis) for three costs of adaptation (columns). Level A0 indicates situations in which the criteria for emergence are met right from the initial conditions. Line 2. An in line 1 when the resistancaffects only infection efficiency (IE scenario). Line 3. Shape of the ℛ 0 function. Three shapes are classified: (i) "1S-2S" when the evolutionary attractor * is closer to S than to R (i.e., * < 0), (ii) "1R" when * is closer to R and ℛ 0 is unimodal, and (iii) "2R" when * is closer to R and ℛ 0 is bimodal. Point A0 indicates situations in which the criteria for emergence are met right from the initial conditions. For all panels, all the other parameters are set to their reference values (Table 2) consequences of deploying resistant cultivars are not necessarily aligned, and depend on the pathogenicity trait targeted by the plant resistance genes. From an evolutionary perspective, the emergence time of the strategy best adapted to the R cultivar tends to be shorter when the resistance affects infection efficiency (IE scenario) than when it affects sporulation (SP scenario). In both cases, the emergence time is maximal for an intermediate proportion of R cultivars, as previously reported by Papaïx et al. (2018). In contrast, from an epidemiological perspective, epidemiological control is always greater in the IE scenario than in the SP scenario, as already observed by Iacono et al. (2012); Rimbaud et al. (2018). These general rules are common to all theoretical studies over the last decade investigating the epidemiological and evolutionary effects of deploying quantitative resistances in agro-ecosystems, supporting the robustness of these conclusions.

| Bridging the gap between modeling approaches
All these previous studies account for the interplay between epidemiological and evolutionary dynamics, but they are based on different modeling frameworks, hypotheses, and parametrizations. One of the strong points of our general modeling framework is that it allows us to bridge the gap between different modeling traditions. First, we can investigate both the short-and long-term epidemiological and evolutionary dynamics of the host-pathogen interaction (Day & Proulx, 2004). The short-term dynamics are investigated numerically, but the long-term analysis is analytically tractable and can predict the outcome of pathogen evolution. In contrast to most studies in evolutionary epidemiology, our results are neither restricted to rare mutations, as in the classical adaptive dynamics approach, nor restricted to a Gaussian mutation kernel (see also Mirrahimi (2017)).
When the trait distribution is sufficiently narrow, we expect the population to concentrate around the attractors predicted by adaptive dynamics. However, the full model is required to quantify the speed of evolution, and, notably, the time to the emergence of an adapted pathogen strain, which is often as important for practical purposes than the characterization of the potential long-term outcome. We also note that "fat-tailed" kernels allowing the inclusion of long-distance dispersal events in the phenotype space could also be considered (Appendix S2).
Second, our analysis makes it possible to consider multimodal fitness functions, which are reminiscent of fitness landscapes in population genetics, and to characterize evolutionary attractors at F I G U R E 5 Epidemiological outcomes. Epidemiological control is estimated by calculating the gain in relative healthy area duration (rHAD), a proxy for crop yield. Line 1. The resistance affects only total spore production (SP scenario). Epidemiological control is characterized as a function of the proportion of the R cultivar at planting (x-axis) and of the relative effectiveness of the resistant cultivar (y-axis) for three costs of adaptation (columns). Line 2. As in line 1 when the resistance affects only infection efficiency (IE scenario). For all panels, all the other parameters are set to their reference values (Table 2) equilibrium through a detailed description of their shapes (number of modes, steepness, and any higher moments with even order). Our results show that these features strongly affect transient dynamics, which, in turn, shape evolutionary and epidemiological outcomes.
In particular, multimodal fitness functions with close local maxima result in the longest time to the emergence of adapted pathogen strains in the SP scenario. Mathematically, they are characterized by a slow convergence to equilibrium (Burie et al., 2020). In practice, it is possible to take advantage of these properties by wisely choosing the proportion of each cultivar in the agricultural landscape. Indeed, the shape of the fitness function in this scenario is basically the sum of the fitness in each hosts weighted by their proportions in the environment (equation 3).
Third, we propose a simple framework based on Gaussian distributions to describe the pathogen phenotypic landscape associated with quantitative resistance. Variation of the relative effectiveness of resistance and cost of adaptation allows us to describe a whole range of quantitative resistance effects, making it possible to consider the continuum between the two main models used for qualitative plant-parasite interactions: (i) the gene-for-gene model in which a parasite strain may have universal infectivity (R eff = 1 and C adp = 0) and i) the matching-allele model in ich universal infectivity is impossible (R eff = 1 and C adp = 1) (Thrall et al., 2016). By exploring the continuum between these extremes (while noting that R eff and C adp can only tend to 1 in our framework), we can investigate in more detail the impact of deploying quantitative resistances on epidemiological and evolutionary outcomes.

| Quantitative resistances as a major driver of pathogen evolution
The possibility that deploying R cultivars in agricultural landscapes will lead to the diversification of an initially monomorphic population and to the long-term coexistence of different pathogen strategies has been investigated with adaptive dynamics by  and by Gandon (2004) in the general context of the evolution of virulence in multihost parasites. These studies highlighted the strong dependence of evolutionary outcomes on the shape of the trade-off curve for pathogen transmission on sympatric hosts. Concave trade-off curves lead to monomorphic evolutionary endpoints, whereas convex or sigmoidal trade-off curves can lead to evolutionary branching. These results can also explain the existence of sibling fungal pathogens (Gudelj, Fitt, et al. (2004)). In this approach, the transmission rate aggregates infection, spore production, and dispersal into a single proxy trait affecting secondary infection from infected to healthy hosts. In contrast, our framework allows us to disentangle the specific effect of the pathogenicity traits that follow one another during infection. Consistent with , we show that dimorphism is possible with a convex or sigmoidal trade-off for host infection efficiencies (IE scenario), whereas the population always remains monomorphic when the trade-off curve is concave. However, our analysis shows that predictions crucially depend on the step of the parasite life cycle affected by quantitative resistance in the R cultivar. The population always remains monomorphic when resistance affects only the sporulation curve (total spore production but also any other pathogenicity traits, such as the duration of the latency period) irrespective of the underlying trade-off curve shape (SP scenario). This is because the SP scenario allows an optimization principle based on ℛ 0 and potential evolutionary attractors are located at the peaks of ℛ 0 (Lion & Metz, 2018).
The evolution of plant pathogens toward generalism or specialism following the deployment of an R cultivar is also a major concern (Croll & McDonald, 2017;Papaïx et al., 2013Papaïx et al., , 2014. We propose a simple binary criterion for classifying the evolutionary attractors as generalist or specialist by comparing their ℛ 0 on the R and S cultivars to a preset threshold ( Figure S4). This highlights not only that high costs of adaptation are a major factor leading to the selection of specialists but also that emergence times and epidemiological control are weakly structured by this classification ( Figure   S4). It remains unclear whether it is better to favor deployment strategies leading to the selection of generalist or specialist pathogen strains, and this issue merits further investigation. The answer depends partly on the preexisting level of pathogen diversification (Papaïx, 2011). However, addressing this question in the context of quantitative resistance will first require the development of new indices measuring the degree of generalism within a specialist/generalist continuum.

| 0 in heterogeneous host environments sharing a common pool of propagules
Many practical epidemiological studies are based on the concept of the basic reproduction ratio ℛ 0 . For instance, van den Bosch et al.
(2008) calculated ℛ 0 for lesion-forming foliar pathogens in a setting with two cultivars but with no effect of the time since infection on sporulation and disease-induced mortality. ℛ 0 is typically calculated with the spectral radius of the next-generation operator (Diekmann et al., 1990), but we follow here a methodology based on the generation evolution operator (Inaba, 2012) to derive an expression for the basic reproduction number ℛ 0 in heterogeneous host populations with any number of cultivars. This expression captures the time scales inherent to the life cycle of plant fungal pathogens with, notably, sporulation varying over time. Capturing such patterns is a challenge in the modeling of plant diseases, as reported by Cunniffe et al. (2015). Moreover, with a common pool of well-mixed airborne pathogen propagules, the function ℛ 0 (x) is an exact fitness proxy for competing strains with potentially different sporulation curves (including the latent period, total spore production, and shape parameters). This is important because the computation of ℛ 0 (x) usually assumes that the invading pathogen enters an uninfected host population. However, more generally, a clear distinction between pathogen invasion fitness ℛ (x, y) and epidemiological ℛ 0 (x) is required for correct discussions of the adaptive evolution of pathogens (Lion & Metz, 2018). Even with a common pool of spores, the optimization principle of ℛ 0 (x) does not hold when infection efficiencies differ between host classes.

| Notes on model assumptions
The model assumes an infinitely large pathogen population.
Demographic stochasticity is, thus, ignored despite its potential impact on evolutionary dynamics (e.g., lower probabilities of emergence and fixation of beneficial mutations, and reduction of standing genetic variation (Kimura, 1962)). In particular, genetic drift is more likely to affect the maintenance of a neutral polymorphism than that of a protected polymorphism when selection favors the coexistence of different genotypes protecting against invasions by mutant strategies (Geritz et al., 1998). The effect of genetic drift depends on the stability properties of the model considered. As our model has a unique globally stable eco-evolutionary equilibrium, genetic drift is likely to have a weaker impact than in models with several locally stable equilibria. Moreover, the large effective population sizes (from 10 3 to 3.10 4 ) reported at field scale for several species of winddispersed, spore-producing plant pathogens (Ali et al., 2016;Walker et al., 2017;Zhan et al., 2001) suggest a weak effect of genetic drift.
The model assumes that the aggressiveness components are mutually independent. However, correlations between traits have sometimes been reported. For instance, Pariaud et al. (2013) observed a positive correlation between the duration of the latent period and the fecundity of a plant fungus. This relationship describes a phenotypic trade-off because a short latent period and a high sporulation probability represent fitness advantages. It can be introduced into the SP scenario assuming that the latent period and the total spore production of a pathogen strain x are linked by mathematical functions, such as the quadratic relationship suggested by Pariaud et al. (2013). Alternatively, correlations between pathogen life-history traits can emerge from the covariance matrix of the (multidimensional) mutation kernel (Gandon, 2004). Another feature of the framework that we did not use is the possibility of tracking the evolution of the quantitative traits of the pathogen with Price's equation Day & Proulx, 2004;Iacono et al., 2012). Indeed, differential equations for the mean phenotype and phenotypic variance of any trait of interest can be derived from model (2).
The model also assumes a unique pool of well-mixed propagules.
Spore dispersal therefore disregards the location of healthy and infected hosts. This assumption, which ensures the one-dimensional environmental feedback loop of the model, is more likely when the extent of the field or landscape considered is not overly large with respect to the dispersal function of airborne propagules. Airborne fungal spores often disperse over substantial distances. Mean dispersal distances range from 100 m to 1 km and, in most cases, longdistance dispersal events are frequent . A spatially implicit framework embedded into integro-differential equations was recently used to describe the eco-evolutionary dynamics of a phenotypically structured population subject to mutation, selection, and migration between two habitats (Mirrahimi, 2017). With the same goal, Rimbaud et al. (2018) used a spatially explicit framework with several habitats embedded in an HEIR stochastic model. It would be interesting to draw on these examples and extend our approach to a spatially explicit environment. Indeed, when dispersal decreases with distance, large homogeneous habitats promote diversification, whereas smaller habitats, favoring migration between different patches, hamper diversification (Débarre & Gandon, 2010;Haller et al., 2013;Papaïx et al., 2013Papaïx et al., , 2014.

| Applied significance
This study gave rise to three key applied findings. The first is the ℛ 0 equation in an environment with several cultivars. This equation is based on a gamma sporulation curve, as documented for several plant fungi (van den Bosch et al., 1988;van den Bosch et al., 1988;Kolnaar & Bosch, 2001;Sache et al., 1997;Segarra et al., 2001), and explicitly integrates the pathogenicity traits expressed during the basic steps of infection (infection efficiency, latent period, sporulation dynamics). As these traits can be measured in the laboratory, this expression of ℛ 0 bridges the gap between plantscale and epidemiological studies, and between experimental and theoretical approaches. It helps address the practical difficulty of estimating ℛ 0 for real populations, as pointed out by Gilligan and van den Bosch (2008). For example, it can be used to compare the fitness of a collection of pathogen isolates (e.g., see Montarry et al. (2010) for an application to potato late blight) or to predict their field structure (e.g., see Durand et al. (2017) for an application to Lyme disease).
The second applied finding is that the choice of quantitative resistance genes operated by plant breeders drives pathogen diversification. This effect is mediated by the parasite life cycle targeted by the resistance. Quantitative resistances specifically affecting different stages of pathogen cycles have been identified in several pathosystems (e.g., Azzimonti et al. (2013), Chung et al. (2010), Jorge et al. (2005)). Moreover, understanding the conditions for maintaining pathogen polymorphism is a long-standing question in disease ecology and evolution, with implications for disease management (Vale, 2013). We show here that the pathogen population always remains monomorphic when resistance affects only sporulation, irrespective of the underlying trade-off curve shape. In contrast, pathogen diversification is possible when resistance affects only infection efficiency. This finding highlights the need for a detailed knowledge of the effect of resistance genes on pathogen life cycle for a full understanding of the impact of the deployment of quantitative resistances on the strain structures underlying the patterns of disease incidence. This finding is important, because quantitative resistances are increasingly being used by plant breeders to deal with the rapidity with which major resistance genes break down (Zhan et al., 2015). Interestingly, similar questions concerning the epidemiological and evolutionary consequences of vaccination strategies have been considered.
Quantitative resistance traits targeting pathogen infection rate, latent period, and sporulation rate are analogous to partially effective (imperfect) vaccines with anti-infection, antigrowth, and antitransmission modes of action, respectively (Gandon et al., 2001).
Similarly, the proportion of the R cultivar deployed is analogous to vaccination coverage in the population; in the vaccination context, the relative effectiveness of resistance is termed "vaccine efficacy" and the relative cost of adaptation is termed "cost of escape" .
Finally, using a different theoretical approach, we confirm the finding of Papaïx et al. (2018) that "there is no silver bullet deployment strategy" simultaneously maximizing epidemiological and evolutionary outcomes (Rimbaud et al., 2021). Accordingly, as the different stakeholders involved in plant resistance management may pursue objectives that are not always compatible (e.g., growers, breeders), they should all be involved in defining the best management strategies.

ACK N OWLED G EM ENTS
We thank Ludek Berec, Eva Kisdi, Loup Rimbaud, and Julien Papaïx for comments and suggestions to improve the manuscript.

AUTH O R S' CO NTR I B UTI O N S
RDD, AD, JBB, and FF planned and designed the research. RDD, AD, JBB, and FF developed analytic tools and conducted numerical experiments. SL and QR critically revised the manuscript and contributed to its writing. All the authors gave final approval for publication and agreed to be held accountable for the work described.

DATA AVA I L A B I L I T Y S TAT E M E N T
The MATLAB codes used to simulate the model and generate