Rotating and stacking genes can improve crop resistance durability while potentially selecting highly virulent pathogen strains

Overview of the model

This model simulates the population and evolutionary dynamics of different pathogen strains, as they interact with different crop resistant cultivars planted in a single field over successive years. We assume one cultivar is planted each year and we consider a field divided into a finite number m of spatial units (representing limited spaces for infections, or potential lesion sites), in which the spatial aspect is implied rather than explicitly represented. For each year during the cropping season, a number of pathogen spores are released from the infested crop residues, it then lands on the crop plants leading to infections (Fig. 1). These infections are apportioned between the different pathogen strains depending on their previous abundance and interactions with the crop cultivar. At the end of the year, during the non-cropping season, these strains are assumed to sexually recombine in the crop residue. The number of spores released and the number of infections are considered as random variables. We denote these both quantities with an uppercase letter (for example N) in general sense, while their particular realisation or draw in the simulation will be noted with a lowercase letter (for example n). The model was developed using the R Language and Environment for Statistical Computing⁴⁹.

We applied the general model described above to a specific situation with four genes of interaction where there are four different resistance genes that may or may not be deployed within each crop cultivar, and four virulence genes that may or may not be present within each pathogen strain. We assume that the presence of each virulence gene reduces the fitness of a strain independently. Specifically, for each strain i, we define the fitness of that strain (delta _i = delta ^{n_{vir,i}}), where (n_{vir,i}) is the number of virulence alleles present in strain i, and (delta) is a fixed model parameter with potential values between zero and one (Table 1).

We first set the model parameter values to define a baseline situation where there is a relatively small fitness penalty for virulence alleles (i.e. (delta) is very close to 1, where the value 1 means no penalty); the pathogen has a relatively low ability to reproduce if it does not carry effective virulence genes (i.e. low value for (epsilon), in this baseline situation equal to 0.05); the initial virulence allele frequency (Init.freq) is relatively low, reflecting a low historical selection pressure and lastly the initial quantity of pathogen (Init.path) is also low at (10 %) of carrying capacity (Table 1). We then considered and compared four different strategies for rotating resistant crop cultivars:

1. S1.

   No rotation, the same cultivar with only one gene of resistance is employed every year;

2. S2.

   A cultivar with a single gene of resistance is employed each year, and the gene of resistance in the cultivar is changed every year, giving a 4-year rotation;

3. S3.

   A cultivar with a single gene of resistance is employed each year, and the gene of resistance in the cultivar is changed every 5 years, giving a 20-year rotation; and

4. S4.

   A cultivar with two genes of resistance (i.e. pyramided resistance) is employed each year, and the genes of resistance in the cultivar are changed every year, giving a 2-year rotation

We then investigated how different parameterisations of the model would interact with the selected rotation strategies. We develop four cases in addition to the baseline case described above:

1. Case 1.

   Baseline scenario (Table 1).

2. Case 2.

   Baseline scenario, except for Init.freq which was increased from 0.05 to 0.5.

3. Case 3.

   Baseline scenario, except for Init.freq which was increased from 0.05 to 0.5 and (delta) which was decreased from 0.9 to 0.7.

4. Case 4.

   Baseline scenario, except for Init.freq which was increased from 0.05 to 0.5, (delta) which was decreased from 0.9 to 0.7 and (epsilon) which was increased from 0.05 to 0.5.

5. Case 5.

   Baseline scenario, except for (epsilon) which was increased from 0.05 to 0.5.

Genetics

Cultivar and pathogen strain are both defined through their genotype being restricted to a specific set of interaction genes (loci) related to resistance (for the cultivar) or virulence (for the pathogen). Each gene has two versions (alleles): virulence or avirulence allele for the pathogen and resistance or susceptibility for the cultivar. Virulence and resistance are represented with a 1 and avirulence and susceptibility are represented with a 0 (Fig. 2). If we call ({mathcal {I}}) the set of strains and if (nu) genes of interaction are involved, then the total number of strains will be (left| {mathcal {I}}right| = 2^{nu }). During the infection process, after pathogen spores land on the cultivar, an interaction factor (beta (i,c)) defines the relative rate at which strain i can reproduce within a field of cultivar c, for each strain and cultivar combination (Fig. 2). We consider that a strain overcomes the cultivar genotype when the strain has a virulence allele for every resistance allele of the cultivar (Fig. 2), in which case (beta (i,c)=1), indicating maximum reproduction rate. Otherwise, if the strain does not have a virulence allele for every resistance allele of the cultivar, (beta (i,c)= epsilon), where (epsilon) is a model parameter with constant value (0 le epsilon < 1), indicating a less-than-maximum reproduction rate. As such, (epsilon) is the model parameter modifying the growth and reproduction of pathogen strains not carrying multiple virulence alleles (e.g. 0100) and/or an avirulent pathogen strain (e.g. 0000) (Fig. 2). Accordingly, lower (closer to 0) (epsilon) values represent reduced ability to grow and reproduce in pathogen strains with increasing number of avirulence alleles. Moreover, any strain i with one or more virulence genes is also assumed to suffer a fitness penalty (delta _i) depending on the number of genes involved. Together these interaction factors make a cultivar-strain interaction matrix (B = (beta (i,c))). This code and method for modelling resistance and virulence interactions (without fitness penalty) is similar to those in previous studies^{48,<a data-track="click" data-track-action="reference anchor" data-track-label="link" data-test="citation-ref" aria-label="Reference 50" title="Lo-Pelzer, E. Modelling the effects of cropping systems and their spatial distribution on phoma stem canker on winter oilseed rape and the adaptation of fungal pathogen populations responsible for the disease (Leptosphaeria maculans) to cultivar resistances. Life Sciences, AgroParisTech, English.. (2008).” href=”https://www.nature.com/articles/s41598-020-76788-7#ref-CR50″ id=”ref-link-section-d44078e2012″>50}.

Initial genetic structure of pathogen population

At the start of each case, we define the initial proportion of each pathogen genotype using the equation:

where strains.init is the initial proportion of each pathogen genotype; Init.freq is the frequency of the virulent genes as set by each case and nr is the number of virulent genes present in a given pathogen genotype. We then used a random Poisson distribution generator (rpois function from the stats package in R) to obtain the initial number of spores for each pathogen genotype, where the mean of the Poisson distribution is the proportion of a given pathogen genotype multiplied by the pre-determined pathogen load (Init.path, Table 1).

Model dynamics

The annual dynamics (Fig. 1) can be divided into two main phases: the phase of parasitic activity, representing events occurring through the cropping season, and the phase of dormancy, representing events occurring between the cropping seasons. During the phase of parasitic activity, the pathogen produces spores which are spread both through the air (sexual ascospores) and via water splash (asexual conidia). These spores may then infect leaves and stems of the cultivar, resulting in new lesions of different strains. During the phase of dormancy, the pathogen remains within the infected crop residue and sexual recombination occurs. These processes are modelled with four steps, three for the parasitic phase and one for the dormancy phase.

Total amount of spores released

First, the model generates the amount of pathogen spores of each strain that is released, using the equation:

where (lambda _{{ released},i}(t)) represents the expected dispersed propagule (spore) pressure in the field due to strain (i in {mathcal {I}}) during the year t, the parameter (alpha) represents the normal rate of growth for the pathogen from 1 year to the next, and (n_{{ recombined},i}(t-1)) represents the number of spatial units or locations infected by the strain i at the end of the previous year and after genetic recombination. The actual quantity of pathogen strain i released in the current year, (n_{{ released},i}(t)) is then simulated as a Poisson random variable:

The infective pressure (lambda _{{ infected},i}(t)) is then calculated as:

where (beta (i,c(t))) is the interaction factor between the strain i and the cultivar c(t) i.e. the cultivar grown in year t, and (delta _i) is the fitness penalty for the particular strain i.

Total number of infections

Second, the model calculates the total number of infected sites, following a binomial distribution:

where (rho _i(t)) is the probability that a particular given location (among the m possible locations in the field) during year t, will have a given spore from strain i fall down on it and cause a lesion, and thus (displaystyle 1 – prod nolimits _{i = 1}^{2^{nu }} (1- rho _i(t))^{n_{{ released},i}(t)}) represents the probability that at least one of the (displaystyle n_{ released}(t) = sum nolimits _{i = 1}^{2^{nu }} n_{{ released},i}(t)) spores produces a lesion. This equation can be justified in more detail as follows:

We assume that a spore will fall on any of the m specific locations with the same probability independently of its infection capabilities. The number of locations m is assumed to be the same for all years whatever the cultivar grown and thus, this probability is independent of the time dimension. Next, we assume that the probability that a spore will induce an infection depends on the interaction factor between the crop cultivar genotype and the pathogen strain (beta (i,c(t))) together with the fitness penalty for that strain (delta _i) . These assumptions mean that:

Number of infections for each strain

Third, the number of infections of each strain is derived from the total number of infections depending on genetic interactions between each strain and crop cultivar being employed that year. Specifically, the total number of infections (N_{ infected}(t) = n_{ infected}(t)) is distributed among the different strains using the multinomial distribution:

where (displaystyle lambda _{{ infected}} (t) = sum nolimits _{i = 1}^{2^{nu }} lambda _{{ infected},i}(t)). The number of infected sites due to strain i, without no loss of generalities, follows then the binomial distribution (displaystyle {mathcal {B}}left( n_{ infected}(t),frac{lambda _{{ infected},i}(t)}{lambda _{{ infected}}(t)}right)).

Genetic recombination

The fourth step involves simulating the process of sexual recombination, where new quantities of each strain are generated based on the previous quantities of each strain. At the end of the year t, we calculate the frequencies (f_j(t)) of each virulent version of each gene from the different genotypes of strains in the crop stubble. We let the genotype of any new spore be represented by a random vector (displaystyle G_i(t) = left( G_{i,1}(t), ldots ,G_{i,nu }(t)right)), where each (G_{i,j}(t)) is a Bernoulli random variable (displaystyle {mathcal {B}}(1,f_j(t))). This vector representation of genotype follows the coding illustrated in (Fig. 2). Assuming that strains recombine independently gene by gene, the probability that (G_i(t)) will be a particular genotype (displaystyle g_i(t) = left( g_{i,1}(t), ldots ,g_{i,nu }(t)right)) is given by:

We can also confirm that across all possible genotypes these probabilities sum to one:

If we shorten the notation for (P(G_i(t) = g_i(t))) to be (p_i(t)) then we can define the recombined version of infected numbers of units of each strain with the following multinomial distribution:

Poisson, binomial and multinomial distribution

In plant pathology, it is often relevant to model infections by a random variable. Let’s imagine a released spore flying in the air, we can say that this spore will land on a specific leaf and infect it with a given probability p, then it won’t with probability (1 – p) because these are the only two possible events. We can define Y a random variable to model the situation. If we say the event ({Y = 1}) represents the success of the event (landing and infection) and ({Y = 0}) represents the failure, with this definition we say that Y follows a Bernoulli distribution. The values attributed to the variable depending on the events allow the following generalisation: If we consider n spores, each of them realizing an infection on a specific plant area they fell on with the same probability p, then we can associate to each spore a Bernoulli distribution (Y_i) where (i in {1,ldots ,n }). If we are interested in the total number of infections occurring on this leaf, assuming the fact that they will happen independently of each other, we can model this situation by the variable (displaystyle S = sum nolimits _{i = 1}^n{Y_i}), called binomial variable. We can also denote briefly (S hookrightarrow {mathcal {B}}(n,p)), where n represents the number of events and p the probability of success of each event. Moreover, the Bernoulli variable Y is related to binomial distribution in the way that we can write (Y hookrightarrow {mathcal {B}}(1,p))^51,52.

Usually it is more likely to model such events by a Poisson law rather than binomial law^53,54. When the number of events is so big that we can approximate it by infinity, and when the probability of success of each event is very small, close to zero, it is possible to link both Poisson and binomial distribution through their respective expectations. So if (lim nolimits _{begin{array}{c} n nearrow +infty p searrow 0 end{array}} n*p = lambda ,) then if we define (X hookrightarrow {mathcal {P}}(lambda )), we have (S xrightarrow {text {distribution}} X). Returning to our example, that means that if we have a ‘close to infinity’ number of spores that could fall onto a given plant and infect it with a very small probability p for each of them and still acting independently, we can model the total number of infections by both S or X. Even if there are millions and millions of spores released, this amount is still small compared to infinity, so using X is still a modelling approximation. The use of binomial or Poisson laws depends on the complexity of the situation. For example, if the modeller wants to simulate a model where he anticipates 15 infections, they can use (X hookrightarrow {mathcal {P}}(15)) or (S hookrightarrow {mathcal {B}}(10000,0.0015)) or (S hookrightarrow {mathcal {B}}(1000000,0.000015)).

We consider now a situation where the plant is attacked by a big number of spores, but with different genotypes modifying their ability to infect, some strains being more efficient than others. To model this situation, we can use a vector of variables, each component representing the number of successes due to a specific genotype. We can choose a vector of binomial number or Poisson number. If we consider the case of a threshold in terms of available space to be infected (a maximum number of infections for the plant), such that spores of different strains are competing for those places, we suggest using a vector of random numbers that follows a multinomial law. This distribution derives from the binomial law, although each component is a specified binomial distribution defined from the parameters of the multinomial distribution. But, it is still possible to interpret some of these components via a conditional Poisson distribution.

From binomial to multinomial distribution

The binomial distribution is a particular case of the multinomial distribution. We consider S a binomial distribution of parameters (n, p) counting the number of success of n independent events where the basic probability of success is p. Let U the random variable be defined by (n-S) the number of failures. In the case where S represents the number of infections, U represents the number of uninfected places. The probability to get k infections is given by:

As a result, the probability of having k success is the same that having (n-k) failures. Then the Eq. (9) shows that U follows a binomial distribution with parameters ((n, 1-p)). We can also say that the couple (S, U) follows a multinomial distribution of parameter ((p, 1-p, n)), that we can denote ((S,U) hookrightarrow {mathcal {M}}(p, 1-p, n)). In a more general way, the analogue of the binomial distribution is the multinomial distribution, where each trial results in exactly one of some fixed finite number k possible successes, with probabilities (p_1), …, (p_k) (so that (p_ige 0) for i = 1, …, k and (sum nolimits _{i=1}^k p_i = 1)), and there are n independent trials. Then if the random variables (X_i) indicate the number of times outcome number i is observed over the n trials, the vector (X = (X_1, ldots , X_k)) follows a multinomial distribution with parameters n and p, where (p = (p_1, ldots , p_k)), that we can also write ({mathcal {M}}left( p_1,ldots ,p_n, N = kright))⁵⁵.

From Poisson to multinomial distribution

We consider here a total number of successes (meaning in our example a number of spores that fall on a place and infect it) X being the sum of the infections due to (omega) different strains (X_i) ((1le i le omega)). If we consider that each (X_i) follows a Poisson law of parameter (lambda _i) and that they are all independent, then X follows a Poisson law of parameter (displaystyle lambda = sum nolimits _{i=1}^{omega } lambda _i). The distribution of each (X_i) conditionally to (X = n) is a binomial distribution ({mathcal {B}}(n,frac{lambda _i}{lambda })). We can prove it for all variable (X_j), with (j in {1,ldots ,omega }):

that we obtain using the Bayes formula for conditioning and the use of independence between the (X_i)’s. Then we replace the probabilities by their Poisson values:

Generalizing this result to the random vector of the (displaystyle (X_i)_{1 le ile omega }) for (omega) strains, the distribution of this vector conditionally to the total number X is a multinomial distribution ({mathcal {M}}left( frac{lambda _1}{lambda },ldots ,frac{lambda _n}{lambda }, X = nright))⁵⁵.

Properties of the model

Let (X_1),…, (X_{2^{nu }}) independent random variables such that (X_j hookrightarrow {mathcal {P}}(lambda _{{ infected},j}(t))) for all (j in {1,ldots ,2^{nu }}), we have the following results:

A. When (m rightarrow infty), (displaystyle N_{ infected}(t) hookrightarrow {mathcal {P}}(sum nolimits _{j=1}^{2^{nu }} lambda _{{ infected},j}(t))),

B. For all (j in {1,ldots ,2^{nu }}), (displaystyle N_{ infected, j}(t) xrightarrow {text {distribution}} X_j left| right. sum nolimits _{i = 1}^{2^{nu }} X_i = n),

C. With A and B when (m rightarrow infty), it is equivalent to either simulate (N_{ infected}(t)) then the conditional multinomial vector (displaystyle left( N_{{ infected},1}(t), ldots , N_{{ infected},2^{nu }}(t)right)) conditionally to the realisation (n_{ infected}(t)), or to simulate directly the previously defined variables (X_1),…, (X_{2^{nu }}).

The number of infected sites due to strain j, without any loss of generalities, follows the binomial distribution (displaystyle {mathcal {B}}left( n_{ infected}(t),frac{lambda _{{ infected},j}(t)}{lambda _{{ infected}}(t)}right)). It is important to notice that it is the same law as a Poisson variable with parameter (lambda _{{ infected},j}(t)) conditionally to the realisation (n_{ infected}(t)) of a Poisson variable with parameter (lambda _{{ infected}}(t)). Referring to formula (10), we can see that when the number of sites available for infection goes towards infinity, meaning that (N_{{ infected}}(t)) behaves like a Poisson law of parameter (sum nolimits _{i = 1}^{2^{nu }} lambda _{{ infected},i}(t)), then the variables (displaystyle left( N_{{ infected},i}(t)right) _{1 le i le 2^{nu }}) behave like independent Poisson law of respective rates (displaystyle left( lambda _{{ infected},i}(t)right) _{1 le i le 2^{nu }}).

Proof of the properties of the model

A. With the help of the reminder, we just have to prove this result:

which could be obtained with the mean value theorem⁵⁶. It means that if the total number of places available for infections was unlimited, these infections could be regarded as being Poisson distributed, with infection pressure as defined previously.

We consider the notation of (5), and to simplify the formula we will note: (rho _i = frac{1}{m} . beta _i) and because the result (10) does not depend on time we reduce the notation such that (10) is equivalent to:

and then we want to prove that:

Replacing m by (frac{1}{x}), with (xne 0), the latest equation equals:

We define (displaystyle f_{beta , n_{released}}(x) = prod _{i=1}^{2^{nu }} f_{i,({beta , n_{released}})}(x) = prod _{i=1}^{2^{nu }} (1 – xbeta _i)^{n_{{ released},i}}). Taking into account the fact that

we apply the mean value theorem (⁵⁶) to the derivable function (f_{beta , n_{released}}), we got the following result:

that finishes the proof of point A.

B. The result is immediate knowing the upper reminder concerning the Poisson–Multinomial laws relationship. We just have to take the value of (omega = 2^{nu }).

C. When m is close to infinity, (N_{ infected}(t)) follows a Poisson distribution whose parameter (expectation) is a sum of parameters. A property of Poisson distribution is that the law of a sum equals in distribution the sum of independent Poisson variables with the respective terms. So that we can rewrite B: For all (j in {1,ldots ,2^{nu }}), (displaystyle N_{ infected, j}(t) xrightarrow {text {distribution}} X_j left| right. N_{ infected}(t) = n).

Source: Ecology - nature.com