Multi‐seasonal modelling of plant‐nematode interactions reveals efficient plant resistance deployment strategies

Abstract Root‐knot nematodes, Meloidogyne spp., are soil‐borne polyphagous pests with major impact on crop yield worldwide. Resistant crops efficiently control avirulent root‐knot nematodes, but favour the emergence of virulent forms. Since virulence is associated with fitness costs, susceptible crops counter‐select virulent root‐knot nematodes. In this study, we identify optimal rotation strategies between susceptible and resistant crops to control root‐knot nematodes and maximize crop yield. We developed an epidemiological model describing the within‐season dynamics of avirulent and virulent root‐knot nematodes on susceptible or resistant plant root‐systems, and their between‐season survival. The model was fitted to experimental data and used to predict yield‐maximizing rotation strategies, with special attention to the impact of epidemic severity and genetic parameters. Crop rotations were found to be efficient under realistic parameter ranges. They were characterized by low ratios of resistant plants and were robust to parameter uncertainty. Rotations provide significant gain over resistant‐only strategies, especially under intermediate fitness costs and severe epidemic contexts. Switching from the current general deployment of resistant crops to custom rotation strategies could not only maintain or increase crop yield, but also preserve the few and valuable R‐genes available.

The following Supporting Information is available for this article:   Methods S1 Computation of the season-to-season basic reproduction numbers R 0 for avirulent and virulent nematodes.
Methods S2 Model fitting to experimental data describing the infection dynamics of susceptible tomato roots by avirulent nematodes during a cropping season.
Methods S3 Sensitivity analysis to assess the parameter impact on the healthy root density (a yield proxy) for a susceptible-only strategy over a 15-season time horizon.
Compartmental diagram representing the plant-nematode interaction model for two successive cropping seasons of resistant (superscript X = R) and susceptible (superscript X = S) plants. Healthy plant roots (H X ) are infected by virulent (subscript v) and avirulent (subscript a) nematodes in the soil (P), before becoming latently infected (E = E a + E v ) and then infectious (I = I a + I v ) feeding sites. Between cropping seasons (shaded area), free living nematodes remain in the soil. Parameters are described in Table 1 (main text). Time horizon Ra tio of re s is ta nt pla nt (%) optimal unconstrained: optimal periodic rotation 100 % S 100 % R range and mean , resistant-only (blue), optimal periodic rotation (yellow) and optimal unconstrained (black). Different unconstrained optimal strategies (yielding the same HRD) were identified, so the ratio is represented by its range (shaded area) and its average value (black solid curve). Default parameter values were used (Table 1 in main text).
Methods S1. Computation of the season-to-season basic reproduction numbers R 0 Season-to-season basic reproduction numbers quantify the multiplication rate of a pathogen from the beginning of a cropping season to the beginning of the next season, when pathogen densities are vanishingly low 1 .
To evaluate the basic reproduction number of avirulent nematodes P a on susceptible plants, we assumed that the dynamics of infected feeding sites E a and I a were fast and that the fraction of virulent offspring δ tended to 0. Under these approximations, the dynamics of free living nematodes P a becomes:Ṗ so that: Healthy roots grow linearly in the absence of nematodes so that at time t during the course of a season H S (t) = H 0 + µxt. Therefore: Taking into account between-season survival ϕ, the season-to-season basic reproduction number of avirulent nematodes on susceptible plants therefore reads: One can proceed in a very same way to compute the basic reproduction number of avirulent nematodes on resistant plants, accounting for the absence of free living nematodes produced during such a season. This yields: which is clearly lower than 1. Finally, since virulent nematodes develop similarly on resistant and susceptible plants, their basic reproduction number is the same on both plants. Applying the method described above and taking into account the presence of fitness costs w β and w r , this leads to: Introducing the effective fitness cost w * = 1 − (1 − w β )(1 − w r ), one gets: Methods S2. Model fitting to experimental data Three parameters could not be set from published data: the infection rate (β ), the conversion factor between root biomass and density of feeding sites (x) and the plant growth scaling factor (k). To estimate these parameters, data obtained by Ehwaeti et al. 2 were used. They describe the infection dynamics of susceptible tomato (cv Moneymaker) roots by avirulent nematodes Meloidogyne incognita, which were initially inoculated in the soil at five controlled densities (P 0 (i) with i = 1, . . . , 5). The nematode density in the roots (N data ) was measured after 42 days and 135 days of cultivation. The relative root biomass (B data ), i.e. the root biomass of an infected plant divided by the root biomass of an uninfected plant, was also measured. Both measures were reported as a mean over five replicates, for each initial nematode density. Only final measures after 135 days of cultivation were used to estimate the parameters (measures after 42 days were used to assess the model validity). Model fitting was performed by finding the x, k, β and y parameter values that minimised the distance between the model output and the data using on a weighted square metric: Parameter y was introduced to take into account the extra root mass corresponding to the root-knot galls caused by nematodes. It was estimated along with the three other parameters. N model (i, x, k, β ) and B model (i, x, k, β ) were computed, based on model (Eqn 1, main text) integrated at time t = 135 days, as follows: with the initial nematode density in the soil P a (0) = P 0 (i), except for H S * which was computed with P a (0) = 0, and with given parameters x, k, β and y (remaining parameter values can be found in Table 1, main text). Division by the means made both N and B terms dimensionless so they could be summed. The computation of the parameter values that minimised J was achieved using the optim function of R with the default Nelder-Mead method.

Methods S3. Sensitivity analysis
We performed a global sensitivity analysis to assess the parameter impact on the healthy root density (HRD) for the 1R + 5S optimal periodic rotation strategy over a 15-season time horizon. The method used is based on factorial design and analysis of variance (ANOVA). The HRD was considered as the output (or observation), the parameters as factors.
We first explored the parameter space by varying all parameters, except the nematode infection success ε X y (a Boolean) and the duration of a cropping season τ. We chose three levels per parameter: the default value found in the literature or estimated according to Supporting Information Methods S2 and ±30% variations, except for P 0 for which larger variations were chosen (Table 1, main text). A full factorial design, defined as all possible combinations of the three parameter levels, corresponds to 3 14 = 4, 782, 969 combinations and would have required the same number of simulations to compute the corresponding HRD. To reduce this number, we implemented a fractional factorial design (a subset of the full design), chosen in order to estimate all parameter main effects and two-way interactions. The fractional factorial design was obtained using the PLANOR R package (https://CRAN.R-project.org/package=planor) and consisted of 2187 parameter combinations, yielding as many simulations.
By means of an ANOVA, we then proceeded with the observed variance decomposition and estimated the sum of squares associated with each factorial term, main effect SS i or the two-way interaction SS i, j . According to the sparsity-of-effects principle, a system is usually dominated by main effects and low order interactions, so neglecting higher interactions can still provide good estimates. Denoting by SS T the total sum of squares, the total sensitivity index of parameter p i is defined as follows: It represents the fraction of the output variability explained by parameter p f . We used the MULTISENSI R package (https://CRAN.R-project.org/package=multisensi) for this analysis. The results of the global sensitivity analysis are shown in Fig. S3. Interactions had little impact, so parameter main effects largely explained the output variability. Parameters that most influenced the HRD were the infection rate β (tSI β = 30% ), the nematode reproduction rate r (tSI r = 28%), as well as the nematode mortality rates in the soil η (tSI η = 23%) and in the plant α (tSI α = 18%). The other parameters had little impact on the output.