Maxime Garnault

Clémentine Duplaix

Pierre Leroux

Gilles Couleaud

Florence Carpentier

Olivier David

Anne-Sophie Walker
email: anne-sophie.walker@inra.fr

SPATIOTEMPORAL DYNAMICS OF FUNGICIDE RESISTANCE IN THE WHEAT PATHOGEN ZYMOSEPTORIA TRITICI IN FRANCE

Keywords: Mycosphaerella graminicola, septoria leaf blotch, spatiotemporal dynamics, fungicide resistance dynamics, DMIs, QoIs, SDHIs

BACKGROUND: Management of pesticide resistance is a major issue in modern agricultural systems, particularly in the context of the broader challenge of reducing pesticide use. However, such management must be adapted to resistance dynamics, which remains challenging to predict due to its dependence on many biological traits of pests, interactions with the environment and pesticide use. We retrospectively studied the evolution of reported resistances to four modes of action (benzimidazoles, quinone outside inhibitors, sterol demethylation inhibitors and succinate dehydrogenase inhibitors); in French populations of the wheat pathogen Zymoseptoria tritici.

RESULTS: We used statistical models to analyse the

Performance trial network dataset (2004-2017; ~70 locations in France yearly). They highlighted contrasting behaviours between phenotypes, such as: (i) stable spatial distributions and colonization front structures over time; (ii) different frequency growth rates at the national scale and between regions.

CONCLUSION:

We provide here a quantitative description of the spatiotemporal patterns of resistance evolution, for fungicides with several modes of action. Moreover, we highlight some unexpected resistance dynamics in France, with major differences between the north and south. This complex pattern of resistance evolution in French populations is consistent with previous descriptions of dynamics at the European scale. These results should make it easier to anticipate the evolution locally and to improve the management of resistances.

INTRODUCTION

In modern agricultural systems, pesticides, or plant protection products (PPPs), make a major contribution to the maintenance of crop productivity, by reducing the impact of pests on yields and quality. However, human health and environment protection issues have driven national and European authorities to adopt policies limiting pesticide registration and use, as in the Ecophyto II plan in France and the 2009/128/CE directive (www.agriculture.gouv.fr). The emergence of resistance is an additional unintended issue of the widespread use of PPPs in agriculture [START_REF] Palumbi | Humans as the world's greatest evolutionary force[END_REF] . Pesticide resistance is the intrinsic, inheritable ability of some pest genotypes to survive pesticide concentrations that kill or inhibit the development of sensitive genotypes of the same species. It results from the adaptive evolution of pest populations under pesticide pressure, leading to the selection and increase in frequency of the least sensitive genotypes [START_REF] Delp | Fungicide resistance: definitions and use of terms[END_REF][START_REF] Barres | Trends and challenges in pesticide resistance detection[END_REF] . This phenomenon may lead to practical resistance, defined as a decrease in the efficacy of a pesticide in the field. In practice, field resistance may lead to an increase in the number of sprays and the dose of pesticide used to maintain the same level of protection, with a greater toxicological and ecotoxicological impact. The development of resistance management strategies is clearly a major challenge, if we are to adhere to a longterm vision of more environmentally friendly agriculture and to avoid the dilemma of "the Tragedy of the Commons" in which multiple individuals acting independently and exclusively in their own interests deplete a limited shared resource, even though this outcome is clearly in no-one's long-term interest [START_REF] Hardin | The tragedy of the commons[END_REF] . For the effective management of pesticide resistance, we therefore need to consider the ability of particular management strategies to limit the selection of resistance to a given pesticide in a given pest, and the factors determining sustainability. Resistance dynamics results from interactions between selection pressures, and biological and agronomic system factors promoting resistance [START_REF] Grimmer | Fungicide resistance risk assessment based on traits associated with the rate of pathogen evolution[END_REF] . A better understanding of these interactions is a prerequisite for the identification of the most suitable strategies for a given situation. Unfortunately, quantitative information about resistance dynamics is often limited, because the collection of such information requires longterm comprehensive monitoring over large territories, funding and stakeholder coordination.

In this context, we focused on the example of fungicide resistance in Zymoseptoria tritici (formerly Septoria tritici; teleomorph, Mycosphaerella graminicola). Z. tritici is an ascomycete that causes Septoria tritici blotch (STB) on winter wheat. During the cropping season, Z. tritici mostly disseminates over short distances by producing pycnidiospores (generated by asexual reproduction and dispersed by splashing on wheat plants), but it can also disseminate over larger distances through the production of ascospores (produced by sexual reproduction, predominantly on crop debris during the short period between crops, but also on the first senescent leaves of wheat during the cropping season 6 ). Z. tritici is a relevant model system for addressing the issues raised above, for a number of reasons. First, it is a pathogen of considerable agronomic relevance. STB is the most damaging disease of wheat in Europe, causing yield losses of up to 50% [7][8][9] . Second, STB is mostly controlled by fungicides in Western Europe, with up to 70% of all fungicide use linked to STB control in the EU 9 . Third, Z. tritici presents biological traits that facilitate its adaptation to selective pressures, particularly those exerted by fungicides: a large population size, considerable genetic diversity, two modes of reproduction and an ability to disperse over long distances 10 . Finally, this pathogen has developed diverse resistance mechanisms, resulting in a range of contrasting adaptation dynamics that can be analysed.

Z. tritici has developed various extents of resistance to fungicides with four of the five authorised modes of action in France: benzimidazoles (e.g. thiophanate-methyl, acting on microtubules), sterol-14αdemethylase inhibitors (DMIs; e.g. cyproconazole, epoxiconazole, metconazole, propiconazole, prothioconazole, tebuconazole), respiration complex III or cytochrome b inhibitors (QoIs; e.g. azoxystrobin, pyraclostrobin, trifloxystrobin) and respiration complex II or succinate dehydrogenase inhibitors (carboxamides or SDHIs; e.g. benzovindiflupyr, bixafen, boscalid, fluopyram, fluxapyroxad, isopyrazam). The last group of fungicides, for which resistance has not yet developed, is multisite fungicides (e.g. chlorothalonil, folpet and mancozeb), which are known to be less prone to resistance development than molecules acting at single sites 11 .

Resistance was selected and became generalized in the late 1980s for antimicrotubule agents (benzimidazole-resistant phenotype: BenR) and in the early 2000s for QoIs (strobilurin-resistant phenotype: StrR), resulting in field resistance and the progressive abandonment of these modes of action for treating STB in cereal crops in Western Europe. Resistance to SDHIs (carboxamide-resistant phenotypes: CarR) is just emerging in Europe and its frequency remains at or close to zero in most countries, including France. The BenR, StrR and CarR phenotypes are associated with unique point mutations of the genes encoding their targets (i.e. tub2 E198A 12 ; cytb G143A 13 ; many single changes in SdhB or SdhC), with a large effect on phenotype 14,15 . They are associated with qualitative resistance (a disruptive shift or single-step resistance), resulting in a bimodal distribution of phenotypes (individuals either sensitive or resistant), with the resistant phenotype rapidly invading the population 16 . By contrast, DMIs have been used intensively for STB control since the late 1970s; the first resistant strains were selected in the 1980s, and field efficacy has since gradually declined, but this decline has differed between molecules and between sites 17,18 . A continuum of decreased sensitivity phenotypes (triazole-resistant phenotypes: TriR) has arisen through the progressive accumulation of mutations in the cyp51 target gene, with an overexpression of cyp51 due to insertions in its promoter region, with or without enhanced efflux correlated with three possible insertions in the promoter region of the transporter MFS1 [19][20][21][22][23] . This last mechanism induces a multidrug-resistant (MDR) pattern, with cross resistance between DMIs, QoIs and SDHIs 19 . The variable assortment of these independent mutations, possibly favoured by sexual reproduction, which occurs annually in this pathogen, and their selection after three decades of continuous DMI use, has led to highly diverse phenotypes in the field. This diversity accounts for the observed quantitative (progressive or slow-shifting), resistance with a multimodal distribution of phenotypes, succeeding and replacing each other, generally with an increase in resistance factors over time. More detailed information about resistance to the fungicides used to treat STB is provided in Supplementary Information 1.

Several phenotypic or molecular methods for detecting and quantifying fungicide resistance in Z. tritici have been reported and can be used to investigate resistance dynamics. Several studies have described the variation over space and/or time of the mean EC50 (half maximal effective concentration: dose at which 50% of the maximal growth inhibition is observed) and its derivatives ("resistance intensity-based approaches" 16,18,[24][25][26][27][28] ), whereas other studies have described variations in the prevalence of several mutations or phenotypes ("resistance frequency-based approaches" 17,18,24,27,[29][30][31][32][33] ). Resistance intensity-based approaches characterize resistance phenotypes clearly, but require strain isolation, which is labourintensive and time-consuming. These approaches are therefore subject to limitations due to the assessment of relatively small numbers of isolates sampled from the population, whereas other tools (e.g. molecular-based quantitative tools or microbiological tests using spore bulks) can quantify resistance frequency without the need for strain isolation, resulting in a more representative sample. Moreover, many of these studies have reported snapshot analyses of resistance status or changes in resistance over a short period 17,24,26,30,33 (~1-3 years) or in a limited territory 25,[27][28][29]32,33 (~2-100 locations), using about 10-50 samples per date or location. Thus, in the studies published to date, the principal limitations to the description of resistance dynamics have been sample size and spatiotemporal scale, whereas the sample selected should be as exhaustive as possible, to make it possible to obtain quantitative descriptions of resistance evolution.

Mathematical modelling can also be used to describe and to predict resistance, as it allows to compare the relative efficacy of different resistance management strategies. Most of the models reported to date describe the selection phase once the resistance has already emerged 34 , but provide theoretical predictions of its evolution, and only one of the published models took field data into account 35 . Far fewer studies have focused on resistance emergence [START_REF] Grimmer | Fungicide resistance risk assessment based on traits associated with the rate of pathogen evolution[END_REF]36,37 , but the first detection of resistance is now predicted reasonably well. Finally, these models focused on the temporal dimension and rarely take spatial heterogeneity into account [38][39][40][41] , even though such heterogeneity of resistance is frequently observed in Z. tritici populations, especially at the scale of regional/national surveys 17,18,26,28,31,32,42 . This may result from the regional heterogeneity in crop environment but also in pesticide use. There is, therefore, a need for comprehensive biological datasets over appropriate and heterogeneous spatiotemporal scales and for models specifically designed for the quantitative description and prediction of resistance evolution in heterogeneous environments.

We retrospectively investigated the evolution of fungicide resistance in French populations of Z. tritici. We used a "frequency-based approach" comparing resistance dynamics between phenotypes for various fungicides corresponding to the different modes of action used for STB control, over wide temporal and regional spatial scales. We used the dataset provided by the "Performance" network, which has operated in France since 2004 providing the annual resistance status in the populations of Z. tritici sampled on winter wheat. A detailed statistical analysis was performed to study how phenotypes resistant to benzimidazoles, QoIs or DMIs dispersed spatially, the speed of such dispersion and differences between regions.

MATERIALS AND METHODS

POPULATION SAMPLING

PERFORMANCE NETWORK DATASET

Most of the Z. tritici populations used for the analysis were provided by the "Performance" network. This trial network includes per year 37-70 very diverse partners, including technical institutes, advisers, merchants, farm cooperatives and chemical companies, and is run by the ARVALIS-Institut du Végétal (technical institute) and INRA BIOGER (for the resistance analyses). According to a common protocol applied annually between 2004 and 2017, this network recorded the efficacy and yield associated with 4-10 sets of fungicide spraying conditions, which were compared with unsprayed plots, in winter wheat. Trials consisted in a randomized block design with 3-4 replicates and a minimum of 20 m 2 for elementary plots. Then, the frequency of resistant phenotypes (see Section 2.2) was determined by analysing infected leaves collected from the plots. Administrative information (e.g. supplier, site and sampling date) and agronomic details (e.g. cultivar, leaf layer sampled) were also recorded. The network conducted 1029 trials in total between 2004 and 2017, throughout France. Sampling matched the pattern of winter wheat cultivation (Figure 1). We used resistance data of control plots only for this study. Indeed, as these plots were not treated with fungicides, their pathogen populations (i.e. the frequencies of the different resistances) can be considered representative of those in the surrounding environment. From 2006 to 2011, control plots were sampled twice: the first sample (T0; n=332 trials) was collected just before the first spray (at about stage Z32 or stem elongation phase, in mid-April to mid-May) and the second (TNT; n=719 trials) 20-30 days after the last spray (at about stage Z39-Z55 or booting and heading stages, in mid-May to mid-June). From 2004 to 2005 and from 2012 to 2017, only the later of these two samples was collected.

ADDITIONAL DATASETS

Some resistant phenotypes arose before 2004. We therefore retrieved, for this analysis, data from previous INRA monitoring programs and from the AFPP network (network of the French Association for Plant Protection). Additional information on resistance frequency was then considered for 1980, 1986, 1991, 1994, 1997 and 1999. For these years, mean resistance frequency was calculated from the observation of 10-20 plots annually, distributed throughout the area under wheat. These data were used to plot changes in resistance status over time, together with the Performance dataset, but were not considered in statistical models (see section 2.4).

RESISTANCE CHARACTERISATION

FUNGICIDES, MEDIA AND CHEMICALS

Fungicides were used as technical-grade compounds at one or two discriminating concentrations, and were kindly provided by their manufacturers: boscalid (1 ppm), carbendazim (5 ppm) and epoxiconazole (0.2 ppm) from BASF Agro (Germany); bixafen (0.5 ppm), prochloraz (0.005 ppm) and tebuconazole (0.05 ppm) from Bayer CropScience (Germany); azoxystrobin (0.5 ppm) and pyrifenox (0.01 and 0.15 ppm) from Syngenta Agro (Switzerland); triflumizole (0.032 ppm and 0.4 ppm) from Nippon Soda (Japan); terbinafine (0.01 ppm) from Sandoz (Switzerland). Tolnaftate (0.25 ppm) was purchased from Sigma Aldrich, St Louis (USA). These fungicides, at the doses indicated, formed a set of discriminatory doses, which were used to test the population. These doses were chosen on the basis of the dose-response curves for individual strains, in which the resistance-associated mutations of target genes were characterised. They were validated by showing that only strains bearing the resistant allele could grow in their presence. Inhibitors of sterol biosynthesis were selected so as to reveal the negative or positive crossresistance specific to some genotypes.

Compounds were dissolved in 80% ethanol and added to molten medium after autoclaving. The concentration of ethanol in the medium, including that used for the controls without fungicides, was 0.5 mL L -1 . Penicillin (5 ppm) and streptomycin (5 ppm), both purchased from Sigma Aldrich (St Louis, USA), were added to the medium for all conditions tested, to prevent contamination with bacteria present on infected leaves.

The PG nutrient medium contained 10 g L -1 glucose, 2 g L -1 K2HPO4, 2 g L -1 KH2PO4 and 12.5 g L -1 agar (pH = 6.3). The medium, to which the antibiotics and fungicides were added at their discriminatory doses, was dispensed into 5.5 cm diameter Petri dishes.

MEASUREMENT OF RESISTANCE FREQUENCY

A "sample" was defined as a population constituted from pycnidiospores produced by 30-40 infected leaves, corresponding to the 3-4 repetitions of the control and kept dried at 4°C until use, one to three months after its collection. Leaves were cut into 2-3 cm pieces with scissors, soaked in 10 ml sterile water for ten minutes and then shaken vigorously. The leaf fragments were then removed. The concentration of pycnidiospores in each bulk sample was estimated with an haematocytometer and adjusted to 5 x 10 6 spores ml -1 . For each population, 300 µl of the bulk spore suspension was spread over a set of Petri dishes, each dish containing a distinct discriminatory dose of fungicide, previously optimized on characterized pure strains to allow the growth of a given resistant genotype, or group of resistant genotypes (see below and Supplementary Information 2). Petri dishes were then incubated in the dark at 17°C. After 48 h of incubation, the dishes were observed under a light microscope (magnification x200). For each discriminatory concentration, sensitive pycnidiospores did not germinate or produced only short germ tubes, whereas resistant pycnidiospores produced long germ tubes (>40% control). Taking into account the percentage of the spores germinating in the control, we estimated the proportion of spores resistant to a given fungicide in the bulk with manual counts, based on ~100 spores in each Petri dish. We calculated the proportions of several phenotypes in the spore mixture. This protocol was used in every experiment during 2004-2017, which ensures the standardization and uniformity of the data. A full description of these phenotypes and associated genotypes is available from Leroux & Walker 19 , but, in short, this test was found to be able to distinguish the following phenotypes, or groups of phenotypes:

-BenR phenotype (specific resistance to antimicrotubule agents): growth in 5 ppm carbendazim -StrR phenotype (specific resistance to QoIs): growth in 0. As correlations are established between each phenotype (defined as the ability to grow under a set of discriminatory concentrations) and a genotype or group of genotypes (Supplementary Information 2) and yearly updated with new genotypes, these phenotype groups can be used as proxies, to follow the evolution of specific mutations of target genes associated with resistance. Indeed, some mutations, or combinations of mutations, were specific to our phenotype groups.

STRAIN ISOLATION AND GENOTYPING

Each year, 50-100 strains were isolated from randomly selected populations (at about 20 sites). Infected leaves were placed in a damp chamber for 24 h to promote cirrhus exudation. Cirrhi were collected with a sterile needle and subcultured on malt-yeast agar medium (MYA: 20 g l -1 malt extract, 5 g l -1 yeast extract and 12.5 g l -1 agar; pH = 5.9) three times, until purity was achieved. Isolated colonies with long germ tubes were also collected from media supplemented with the discriminatory concentrations described above, with a sterile needle, under a stereomicroscope. Colonies were then purified as previously described. Pure strains were stored as spore suspensions in 20% glycerol at -80°C.

DNA was extracted and cyp51 sequences were established as described by Leroux

ESTIMATION OF FUNGICIDE USE

STATISTICAL METHODS FOR ASSESSING RESISTANCE DYNAMICS

GENERAL FEATURES OF STATISTICAL MODELS

The frequency of each resistant phenotype in Z. tritici populations was considered as the dependent variable in our models. We excluded only the TriR1-TriR3 and CarR phenotypes, for which we had too few non- 1). We used a logistic regression framework where the observations followed a binomial distribution of probability p and where the logit of p

Phenotype

First year observed

Number of years observed

Number of regions monitored

Number of regions monitored yearly

Number of observations

(𝑙𝑜𝑔𝑖𝑡(𝑝) = ln ( 𝑝 1-𝑝
)) is a linear combination of explanatory factors. Three models were studied here: a SPATIAL MODEL to determine the spatial structure of each resistant phenotype; an ANOVA MODEL to evaluate the variability of frequencies among years and regions; and a DYNAMIC MODEL to model the dynamics of resistant proportions. The SPATIAL MODEL was considered in a frequentist framework, whereas the ANOVA and DYNAMIC MODELS were analysed in a Bayesian framework.

We considered three main factors in these models: year, region and sampling date. For sampling date, TNT (early summer samplings) was used as the reference (i.e. this parameter was set to 0), because TNT was observed every year in the dataset. We therefore actually estimated the difference between TNT and T0 (mid-spring samplings) effects within a year.

In this study we focused on periods and locations for which resistance was evolving (i.e. after emergence and before fixation). Thus, we selected appropriate years and regions among the dataset. For each phenotype, we defined a time range gathering years for which at least a quarter of the regions were monitored with frequencies differing from 0% and 100% (see time ranges in Table 1). In ALS, FCO, LAR and LIM regions, only one trial per year was conducted over three years or less, unlike in other regions where three trials were conducted over at least five years. For the ANOVA and DYNAMIC MODELS, these four regions were removed from the dataset for all the phenotypes. Then, for each phenotype, we selected regions whose frequency data differed from 0% and 100% during at least a quarter of the years of the previously defined time range (see the number of remaining regions in Table 1). Finally, this selection only removed a bit less than 5% of the observed data that were different from 0% and 100%.

Due to the logit link function involved in these models, values of 0 and 100 were overweighted, making parameters difficult to estimate for several phenotypes (e.g. when resistance was emerging or when resistance was completely generalised). We therefore added 0 and 100 inflation parameters 44,45 to the last two models. These two inflation parameters are linked to the proportions of 0 and 100 values observed in the dataset.

The models described below were run for each phenotype, with over-dispersion parameters in the binomial distribution for the ANOVA and DYNAMIC MODELS. Statistical analyses were performed with R software 46 and the rjags package 47 for Bayesian statistics. Parameter estimates were considered significant at the 5% level (or, alternatively, the 2.5% or 0.1% level) if 0 lay outside the 95% (97.5% or 99.9%, respectively) confidence interval for the parameter. Qualitative factors were tested at the same significance levels (5%, 2.5% and 0.1%), using Wald tests 48 . The goodness-of-fit of models was assessed using the coefficient of determination (R²) which varies between 0 and 1. We used the following definition of the R²:

𝑆𝑆𝑀 𝑆𝑆𝑇
, where SSM corresponds to the explained sum of squares (squared differences between predictions and the average of observations), and SST corresponds to the total sum of squares (squared differences between observations and the average of observations). The closer its value is to 1, the better the goodness-of-fit will be.

MODELLING THE SPATIAL STRUCTURE OF RESISTANCE

This model was designed to visualize the spatial clustering of resistances over time on individual maps. We used a geostatistics method to interpolate from sparse sample data by kriging to produce maps along a regular grid 49 . This method principle is that « data that are close together are usually more correlated than those that are far apart » 50 . The spatial dependency of variables is modeled by a variogram (a function only dependent on distance). We estimated variogram parameters using a Matérn space covariance matrix 51 modelling the covariance among observations in a generalized linear mixed model 52 taking year and sampling date into account as qualitative factors. Using the estimated variogram, we obtained by kriging an unbiased linear prediction of minimum variance 50 of resistance frequencies over the monitored area for a given year. We then used the SpODT method (spatial oblique decision tree) to estimate lines that delimit areas as different as possible 53 . Based on a non-parametric method of classification, this algorithm partitions the area by a line that gives the two spatial classes which maximize the interclass variance. This analysis involved the use of two packages: spaMM 54 for spatial interpolation, and SPODT 53 for spatial partitioning, and is referred to hereafter as the "SPATIAL MODEL".

MODELLING VARIATIONS OF RESISTANCE FREQUENCY BETWEEN YEARS AND REGIONS

This second model was designed to identify specific years or regions deviating significantly from an estimated national mean resistance frequency. This model ranked years and regions according to their observed resistance frequencies. Year and region were considered as qualitative factors. The interaction between year and region, and the impact of sampling date in control plots (T0 vs. TNT) were also included in this model. The statistical form of the model is given in Supplementary Information 3.1 and this model is referred to hereafter as the "ANOVA MODEL". For each phenotype, the estimated model parameters permit to compute adjusted mean i.e. to compensate for data imbalances. Thus, the interannual and national mean frequency is given by (1-π0-π100)*(logit -1 (µ) + π100, where logit -1 (x)= exp(x)/(1+exp(x)) and µ is the intercept of the model. The mean frequency for T0 (mid-spring samplings) is given by (1-π0-π100)*logit -1 (µ-αTNT-T0)) + π100, where αTNT-T0 is the effect of the sampling date. The adjusted frequency for the year j is given by (1-π0-π100)*(logit -1 (µ+αj) + π100 where αj is the effect of the year j. The adjusted frequency for the region k is given by (1-π0-π100)*(logit -1 (µ+βk) + π100 where βk is the effect of the region k. More generally, the frequency estimated for a year j and a region k is (1-π0-π100)*(logit -1 (µ+αj+βk) + π100.

MODELLING THE RATE OF INCREASE IN RESISTANCE

The aim of this third model was to estimate and compare the rates of evolution of the different phenotypes in France and between regions. As above, the model included year, region and sampling date effects, but year was treated here as a quantitative variable. The model was therefore dynamic because the years were ordered. This model provided estimates of growth rates of resistance frequency (hereafter named "growth rate") at national and regional scales. The exponential of these growth rates could be interpreted as the apparent relative fitness of phenotypes or, in other words, how much faster the resistant phenotype evolved than the rest of the population 55 . As we consider a relative growth, estimates will be close to 1 if the frequency of the phenotype remains stable in the population, above 1 if the frequency increases over time and below 1 otherwise. The statistical form of the model is given in Supplementary Information 3.2. This model is hereafter referred to as the "DYNAMIC MODEL".

RESULTS

COMPARISON BETWEEN SAMPLING DATES

For each phenotype, the ANOVA and DYNAMIC MODELS gave similar estimates of the "TNT-T0" parameter. This parameter informs on the intra-annual evolution of resistance frequency in control plots by estimating a constant difference between mid-spring and early summer samples (Table 2 andSupplementary Information 4 and5). Hereafter, we therefore provided the results for the ANOVA MODEL only. Estimates were significantly positive for the StrR, TriR6 and TriMR phenotypes, with increase in frequencies between TNT and T0 of 4.7% (P < 0.001), 8.9% (P < 0.001) and 3.5% (P < 0.025), respectively (Table 2). By contrast, estimates were significantly negative for the TriR2-TriR4 and TriLR phenotypes, with decrease in frequencies between TNT and T0 of -1.87% (P < 0.025) and -2.2% (P < 0.025), respectively (Table 2). Furthermore, a strong positive correlation (Pearson's correlation coefficient of 0.82, P = 0.004) was found between estimated "TNT-T0" parameter and the national growth rate, for all phenotypes. The frequency of resistance between T0 and TNT (intra-annual variation) thus followed the same pattern as general inter-annual variation. The correlation coefficient reached 0.90 (Pearson's correlation coefficient, P = 0.006) if only the phenotypes displaying significant inter-annual evolution were considered (Supplementary Information 5). We focus below on TNT observations, unless otherwise indicated. The table gives the inter-annual national mean frequency and the difference between the adjusted frequency and the mean frequency for each level of each factor. Frequencies are expressed in percentage. Estimates of interactions between region and year are not shown because the factor was not significant for all phenotypes. Empty cells indicate that no data were available for estimation. Factors tests P-values are indicated in bold and italic. P-value thresholds: "•" (P < 0.1), "*" (P < 0.05), "**" (P < 0.025), "***" (P < 0.001).

RESISTANCE TO ANTIMICROTUBULE AGENTS (BENZIMIDAZOLES)

Benzimidazoles were introduced into France in the mid-1970s to control STB (Supplementary Information 1). They have marginal use on wheat with an 𝐻𝐴𝐷 𝐶 index of 0.0081 on average over the last five years (residual use against Fusarium head blight). Over the last decade, the frequency of BenR in Z. tritici populations has remained very high and stable (Figure 2), as already observed in the 1990s 56 . Inter-annual national frequency has been estimated at 90% (Table 2). The stabilisation of this resistance is confirmed by the lack of a significant effect of the year factor (Table 2) and by an estimated increase in frequency not significantly different from 1 at national scale (Figure 3). Thick ticks along the x-axis indicate years for which resistance frequency was assessed by monitoring. Lighter ticks indicate that information was obtained or extrapolated from sources other than the Performance database. BZ: benzimidazoles (antimicrotubule agents); DMI: sterol demethylation inhibitors; QoI: inhibitors or respiration complex III; SDHI: inhibitors of respiration complex II. Phenotypes as described in section 2.2.

A regional effect was detected, with differences between the western and northern regions (P < 0.001; Table 2; Supplementary Information 6.1). Indeed, the estimated frequency of resistance was significantly lower than the inter-annual national mean in the AQU and BRE regions (-10.75% (P < 0.025) and -7.25% (P < 0.025), respectively). Resistance frequency was significantly higher than the inter-annual national mean in the NPC and PIC regions (+4.1% (P < 0.05) and +3.9% (P < 0.025), respectively). The rate of increase in the frequency of BenR was similar to the national trend for all regions, with the exception of PDL, which displayed a greater increase in frequency than the other regions (estimate = 0.07; P < 0.05), and CEN, in which the increase in the frequency of BenR was significantly lower and negative (estimate = -0.07; P < 0 .05) (Figure 3 and Supplementary Information 5). The apparent relative fitness is equal to the exponential of frequency growth rate from the DYNAMIC MODEL, and it represents how much faster the resistant phenotype evolved compared to the rest of the population (see Section 2.4.4).

When the apparent relative fitness is close to 1, resistance frequency tends to be constant in the population. Vertical lines: apparent relative fitness estimates at the national scale for each phenotype; lines are continuous if parameters are significant (P < 0.05), dashed otherwise. Dots: estimates of the regional adjustment of the national apparent relative fitness; s. means significant (P < 0.05), n.s. means non-significant. Regions on the vertical axis are organized according the North-South gradient.

RESISTANCE TO QUINONE OUTSIDE INHIBITORS (QoIs)

QoI fungicides were first used to control STB in France in 1997 (Supplementary Information 1) and the first resistant strains were detected in 2002 56 . QoIs use is currently half as high as in the early 2000s. It decreased substantially between 2005 and 2007 because of their poor efficacy on resistant populations (Figure 2). The remaining use of fungicides with this mode of action is chiefly to control rusts and Fusarium head blight.

The national frequency of StrR increased sharply between 2002 and 2012, as shown by the estimates of the year parameters in the ANOVA MODEL (from -43% in 2004, to + 8.3% (P < 0.001) in 2012), with a frequency of 95% already attained in 2010 (Table 2). The DYNAMIC MODEL confirmed this trend, with a significant and positive growth rate of 0.68 (P < 0.001, Supplementary Information 5), indicating that the StrR population grew twice (e 0.68 ≈ 2) as fast as the rest of the population over the study period (Figure 3). Furthermore, StrR increased in frequency significantly more rapidly than all the other phenotypes (at least P < 0.05, data not shown).

The ANOVA MODEL (Table 2) detected a spatial structure, with northern regions, such as NPC and HNO, having a significantly higher StrR frequency on average over time (+5.2% and +4.3%, respectively; P < 0.05), and southwestern regions, such as AQU and PDL, having a lower StrR frequency on average over time (-13.2% (P < 0.1) and -11.6% (P < 0.05), respectively).

Based on the output of the DYNAMIC MODEL (Figure 3 and Supplementary Information 5), increases in frequency were significantly lower than the national rate for the northern PIC and CHA regions (estimates = -0.35 (P < 0.001) and -0.34 (P < 0.05), respectively), whereas the south-western regions MPY and PCH had significantly higher rates of frequency increase (estimates of 0.53 (P < 0.025) and 0.34 (P < 0.05), respectively). The coefficient of determination (R²) for this model is much higher (0. Resistance frequency followed a north-south gradient between 2004 and 2008, with clear delimitations (advancing front) highlighted by the SPATIAL MODEL (Figure 4 and Supplementary Information 6.2). The speed of the progression line (i.e. the distance between the fronts of two successive years) was 120-145 km/year (depending the timeframe used: 2004-2007 or 2004-2008).

RESISTANCE TO SUCCINATE DEHYDROGENASE INHIBITORS (SDHIs)

RESISTANCE TO STEROL 14α-DEMETHYLATION INHIBITORS (DMIs)

DMIs include imidazoles (e.g. prochloraz), triazoles (e.g. cyproconazole, epoxiconazole, metconazole and propiconazole) and a triazolinethione (prothioconazole, which is metabolised to generate prothioconazoledesthio, a triazole). These compounds have been used to control STB on winter wheat in France since the early 1980s, with new molecules regularly released onto the market (Supplementary Information 1). DMIs are the most widely used fungicides in France compared to the other modes of action. We noticed a steady increase in the use of this mode of action, with the index 𝐻𝐴𝐷 𝐶 ranging from ~2 in the early 1990s to over 3 currently (Figure 2). This observation suggests an increase in DMIs doses and/or in the number of compounds mixed in spraying treatment.

TriLR

TriLR phenotypes include strains with weak specific resistance to DMIs: phenotypes TriR1 to TriR5. These strains were first detected in the late 1980s (Figure 2). As we had too few non-zero frequencies observed for TriR1 and TriR3, we focused our analysis on the TriLR group in general and on phenotypes TriR2, TriR4 (grouped into the TriR2-TriR4 group) and TriR5. The estimated inter-annual national frequency of TriLR according to the ANOVA MODEL was 11% (7.6% for TriR2-TriR4 and 4.5% for TriR5; Table 2). For TriLR, annual estimates decreased from +28% in 2005 (P < 0.001) to -10.3% in 2016 (P < 0.001). This trend was also valid for TriR5 (P < 0.001, for the year factor), but was less pronounced for TriR2-TriR4 (P < 0.1, for the year factor). The national rates of frequency increase from the DYNAMIC MODEL confirmed these trends, as they were significantly negative for TriLR, TriR2-TriR4 and TriR5 (estimates of -0.21, -0.1 and -0.14, respectively; P < 0.001; Figure 3 and Supplementary Information 5).

Spatial differentiation was significant for TriLR and TriR5 (P < 0.05, for the region factor in the ANOVA MODEL), but not for TriR2-TriR4 (Table 2). Similar estimates were obtained for the TriLR and TriR5 phenotypes, with differences between northern regions and central-western regions (Table 2 and Supplementary Information 6.3, 6.5). For example, the region parameters for TriLR were significantly negative for PIC and NPC (-4.6% and -3.9% respectively; P < 0.05). This phenotype was therefore less present in these regions than in other regions over the study period. By contrast, it was significantly more present in the PCH and CEN regions (+6.8% (P < 0.025) and +3.7% (P < 0.05), respectively).

For the TriR2-TriR4 and TriR5 phenotypes, none of the regional rates of increase in frequency differed significantly from the national trend (Supplementary Information 5). Nevertheless, the rate of increase in frequency of the TriLR phenotype was significantly less negative than the national rate in the northern region NPC (estimate = 0.12; P < 0.05), whereas significant counter-selection against this phenotype was observed in the southern region RAL (estimate = -0.13; P < 0.05). Overall, there was counter-selection against TriLR strains nationally during the timeframe of this study, because all regional growth rates (the sum of the national growth rate and regional adjustments, Supplementary Information 5) were negative.

TriMR

TriMR strains display moderate resistance to azoles. They include the TriR6, TriR7 and TriR8 phenotypes, which were selected in the late 1990s (Figure 2). Between 2005 and 2017, the estimated inter-annual national frequency of TriMR was 78% (with 56% for TriR6 and 17.5% for TriR7-TriR8). TriMR, and TriR6 in particular, were the most frequent group of TriR strains between the mid-2000s and 2017 (Table 2).

The ANOVA MODEL provided increasing estimates of annual frequency for TriMR between 2005 (-16.5% 0.83 in 2005; P < 0.025) and 2011 (+7%; P < 0.05). These estimates then decreased to reach -10.4% in 2016 (P < 0.05). This dynamic led us to split the TriMR data into two subgroups (TriMR≤2011 and TriMR≥2011), to estimate growth rates more reliably over these two periods.

The DYNAMIC MODEL estimated a significant positive growth rate of 0.22 for the period before 2011 (P < 0.001), indicating that this population of strains grew 1.25 times (e 0.22 ≈ 1.25) faster than the rest of the population (Figure 3). The rate of increase in frequency was significantly negative after 2011 (estimate = -0.17; P < 0.025).

By contrast, no significant trends were observed for the TriR6 and TriR7-TriR8 phenotypes in the DYNAMIC MODEL (nonsignificant growth rates close to 0, Supplementary Information 5). However, the structure of the TriMR group varied over time. Indeed, the estimates of some year parameters were significant for TriR6, and either positive, when this phenotype was present at higher frequencies than the inter-annual national mean (as in 2009, 2010 and 2011; P < 0.05, about 67% on the frequency scale), or negative, when this phenotype was present at lower frequencies (as in 2013 and 2017; P < 0.05 and P < 0.001 respectively, 44% and 23% on the frequency scale). Conversely, the TriR7-TriR8 phenotypes were present at higher frequencies than the inter-annual national mean in 2013 and 2017 (P < 0.05 and P < 0.001 respectively, 26% and 43% on the frequency scale), and lower frequencies in 2012 (P < 0.025, 9% on the frequency scale).

The SPATIAL MODEL revealed a strong geographical partitioning of TriMR subgroups, with TriR6 strains found mostly in northeastern regions, and TriR7-TriR8 strains found mostly in south-western areas, particularly before 2014 (Figure 4 and Supplementary Information 6.7-8). Furthermore, estimates of TriR6 frequencies (ANOVA MODEL) were significantly greater in the northern regions NPC and IDF, than the inter-annual national mean frequency (+11.4% and +11.2% respectively; P < 0.05). By contrast, estimates for the TriR7-TriR8 phenotypes in the southern regions MPY, AQU and PCH were significantly higher than the inter-annual national mean frequency (+35.3% (P < 0.025), +16.2% (P < 0.05) and +10.2% (P < 0.05), respectively).

Estimates for the northern regions BNO and NPC were significantly below the inter-annual national mean frequency (estimates of -7.4% (P < 0.05) and -7.25% (P < 0.05), respectively). Contrasting selection patterns were observed between regions with the DYNAMIC MODEL (Figure 3 and Supplementary Information 5) for the TriMR phenotypes, particularly for the TriR6 and TriR7-TriR8 subgroups. For example, TriR6 strains were selected, as shown by the positive rates of increase in their frequency, in RAL (estimate = 0.16; P < 0.05) and PCH (estimate = 0.13; P < .05), but counter-selected in MPY (estimate = -0.15; P < 0.05) and CHA (estimate = -0.13; P < 0.025). Similarly, TriR7-TriR8 phenotypes were selected in MPY (estimate = 0.22, P < 0.001) and BOU (estimate = 0.14; P < 0.05), but counter-selected in IDF (estimate = -0.16; P < 0.05).

TriHR

TriHR strains are highly resistant to some DMIs and are associated with the most complex cyp51 genotypes (Supplementary Information 2). They were first detected in our dataset in 2008. Since 2010, their frequency has increased, and the estimated inter-annual national frequency of these strains was 7% (Table 2). According to the ANOVA MODEL, TriHR frequency increased from 1.9% in 2010 (7.1-5.2%, P < 0.025) to 16.4% in 2016 (P < 0.05). The DYNAMIC MODEL confirmed this national tendency, with a significant positive rate of increase in frequency of 0.37 (P < 0.001), indicating that the TriHR population grew almost 1.5 times (e 0.37 ≈ 1.5) faster than the rest of the population (Figure 3).

No significant differences in TriHR frequency between regions were noted over the studied period (Table 2). Similarly, no regional adjustment parameter for the national rate of increase in frequency was significant (Figure 3). Maps from the SPATIAL MODEL are provided in Supplementary Information 6.9.

MULTIDRUG RESISTANCE

In France, MDR phenotypes were first detected in 2008, at about the same time as TriHR strains emerged. MDR is based on an original mechanism of enhanced efflux, giving rise to cross-resistance between DMIs, QoIs and SDHIs. As this mechanism was selected in TriR and StrR backgrounds, resistance factors are high only for DMIs and QoIs, due to the combination of resistance mechanisms. The national frequency of MDR over the study period was estimated at 4.7% (Table 2). No significant year or region effects able to account for the spatiotemporal distribution of this phenotype were identified (Table 2). MDR frequency tended to increase in France, but this trend was not significant (Figure 3 and Supplementary Information 5). Likewise, there was no significant regional effect on the national rate of increase in the frequency of MDR (Figure 3). Maps from the spatial model are provided in Supplementary Information 6.10.

The systematic genotyping of pure MDR strains collected since 2008 revealed that the type I mutation (519 bp insertion in the promoter region of MFS1) was the most frequent allele associated with MDR in French populations (frequency greater than 75% in the dataset; Figure 5).

DISCUSSION

We performed a retrospective study of changes in resistance frequency in Z. tritici populations collected from control plots in the Performance national field trial network. The populations studied were not treated with fungicides, and the resistance frequencies presented here therefore probably underestimate those in farmer's plots treated with chemicals, in which field resistance may occur. We adopted this approach as a means of studying unconstrained resistance dynamics in a population subject to little or no fungicide selection pressure. In the presence of such selection pressure, it would have been difficult to separate out the respective weights of general variation in the population and local selection in treated plots.

Resistance to benzimidazoles (BenR), QoIs (StrR), DMIs (TriR) and SDHIs (CarR) has been monitored and analysed in France since 2004. Our analysis, based on three statistical models, provided a quantitative description of resistance dynamics in time and space.

We found a significant difference in resistance frequencies between the first sampling in spring (T0) and the second sampling in early summer (TNT). The direction of this change within the year depended on the phenotype considered, and was positively correlated with the between-year dynamics of this phenotype. Plots are continually contaminated with extrinsic ascospores during the cropping season 6 , and these ascospores may have contributed to the observed changes in resistance frequency, depending on fungicide selection history within their area of origin. These findings also suggest a possible correlation between fungicide resistance and virulence, as already suggested by Yang et al. 59 and Zhan et al. 60 . Indeed, as recombination is minimal during the cropping season, host selection, leading to an increase in virulent strains 61 may be accompanied by an increase in the frequency of resistance mutations associated with the such genetic backgrounds. These intra-annual variations were taken into account in our models and should be integrated into measurements of selection pressures.

The emergence of resistance in France was consistent with fungicide use, as resistance emerged with a time lag of about six years after the authorisation of QoIs (1996 to 2002) and SDHIs (2006 to 2012). Regarding resistance to DMIs, the first mutant strains (TriLR phenotypes) were observed in the early 1990s, about ten years after registration. TriMR phenotypes, which first contributed to the erosion of DMIs efficacy, were detected in populations since the early 2000s, almost 20 years after the mode of action registration, whilst TriHR and MDR, associated to the greatest resistance factors emerged only in the late 2000s. In our dataset, resistance dynamics differed between phenotypes in terms of the speed at which resistance developed and its spatial progression, often with differences between northern and southern regions.

Resistance

to benzimidazoles (BenR phenotype) is well established in French populations, with a stable inter-annual national mean frequency of about 90%. There seems to be no fitness penalty associated with this resistance, as benzimidazole use is minimal and mostly restricted to Fusarium head blight control. BenR frequency was relatively homogeneous over the entire country, with frequencies slightly lower in the Western regions. By contrast, resistance to SDHIs is just emerging and was detected only at a few sites in northern France, and remains negligible. This slow evolution might be the consequence of the recommendation to use only one SDHI spray per season, massively followed by farmers in France.

Resistance to QoIs (StrR phenotype) emerged in the north of France in 2002, one year after it was first reported in the UK and Ireland, and it has since progressively invaded populations (progression of 120-145 km per year), with a clear delimitation each year between the area invaded by resistant strains and a resistancefree area. This is the first time that such a colonization front structure has been quantitatively described for resistance. QoI resistance, which has reached fixation over the entire country since 2009, was found to be increasing in frequency significantly more rapidly than the other resistances, highlighting the high risk of resistance associated with this mode of action. The resistance dynamics for QoIs can be explained by regular selection pressure (two reduced-dose applications of strobilurin associated to a DMI active ingredient, per year, in the early 2000s, in areas with high disease pressure) and the maternal inheritance of the mutation leading to the G143A change (affecting the mitochondrial gene encoding cytochrome b), enhancing the transmission of this mutation. This resistance remains fixed in French populations (>95%), possibly due to a very small or non-existent cost of resistance and the non-negligible residual use of QoIs against rusts and Fusarium head blight, which may maintain the selection pressure on Z. tritici populations for this phenotype. QoI resistance has been reported to occur in a limited number of independent genetic and/or geographic backgrounds, through at least four recurrent mutations in Z. tritici 31 in Europe. The multiple emergences of such resistance was recently confirmed in North American populations, in which significant mitochondrial genome bottlenecks were observed 62 . We found that the StrR phenotype spread very quickly from northern to eastern regions in France, but that its progression in the south and west was much slower. This anisotropic dissemination highlights a preferred direction of migration, from north-western regions to eastern regions, consistent with earlier results showing that gene flow occurs mostly in a west-to-east direction 31 . The lower levels of southward migration may be explained by several factors, including the prevailing direction of the wind, which carries the ascospores 63,64 ; lower local selective pressure (only one QoI application yearly) due to smaller epidemics in the south; or the occurrence of large durum wheat areas, on which Z. tritici populations have been shown to be genetically different from those on bread wheat populations in analyses based on neutral markers 65 . The propagation speed of 120-145 km/year calculated here is, thus, almost certainly a minimal estimate, because the north-to-south gradient observed in our dataset is not the preferential axis for propagation.

By contrast, DMI resistant strains have been present throughout France since at least 1997 56 , but probably for longer, given that DMI selection pressure began in the late 1970s. In UK, progressive decrease in strains sensitivity has been measured over the last 12 years 16 . Quantitative dynamics analysis revealed that some phenotypes were probably counterselected (negative rates of increase in frequency regularly recorded during the 2004-2017 period, e.g. TriLR, associated with deletions or changes in Cyp51 at codons 459, 460 or 461, or the Y137F or V/C136A variants; TriMR after 2011, associated with the I381V and A379G variants), whereas some others were positively selected (positive rates of increase in frequency regularly recorded, e.g. TriMR before 2011; combined changes in TriHR genotypes, including D134G, V/C136A or S524T; MDR, associated with three mfs1 alleles). This situation is consistent with selective replacement favouring a high diversity of resistant genotypes over the last 30 years. Novel cyp51 resistance genotypes evolved through de novo mutation and intragenic recombination. They are thought to have arisen locally (putatively in UK or Denmark) only once or twice and to have dispersed eastwards across Europe through wind-dispersed ascospores. Regular recombination associated with sexual reproduction, and the continuous use of azoles have led to increases in resistance frequency, particularly for novel genotypes 64 . Similar patterns of diversifying selection and intragenic recombination, leading to the emergence of new genotypes with a selective advantage, have recently been observed in North America 62 . The observed dynamics are consistent with the general gradual erosion of azole efficacy observed in several European countries 16,24,28,29 . Results from a European trial network demonstrated major variations in azole performance across Europe correlated with a clear eastward pattern of decrease in the frequency of all cyp51 mutations, with the exception of the I381V and A379G variants (TriMR genetic background 18 ). Spatial heterogeneity in azole efficacy, correlated with population structure for the frequency of resistant genotypes, was also observed in our dataset at the scale of France. Indeed, from 2005 to 2013, static front structures were detected for TriR6 strains, present mostly in North-East France, and for TriR7 and TriR8 strains, present mostly in the South-West. TriR7 and TriR8 strains were found to be moderately resistant to all DMIs, with the exception of prochloraz 66 , potentially accounting for the higher efficacy of this fungicide in this area (Arvalis-Institut du Végétal, unpublished). This spatial structure may result from a founder effect, bringing TriR7 and TriR8 genotypes from northern Europe to the South of France, where they are adapted to the local environment (e.g. climatic conditions, cultivars, low prochloraz use). This hypothesis should be explored by performing fitness tests and analysing cyp51 diversity between strains from northern and southern populations. In more recent French populations from the North and South, contrasting azole efficacies were also observed, with reduced epoxiconazole, prothioconazole, metconazole and tebuconazole efficacy positively associated with higher frequencies of the CYP51 V/C136A and D134G changes, largely present in TriHR genotypes, in the North, compared to the South, in 2016 18,67 .

The large amount of data made it possible to perform a fine quantitative and exhaustive description of resistance dynamics over time and space in this study. To our knowledge, this is the first population study on Z. tritici including resistance phenotyping to all unisite modes of action used to control STB, related selection pressures, a large temporal scale (13 years) and a fine spatial coverage (including regional scales for the entire wheat growing area in France), in contrast to previously published snapshot analysis. Our resistance frequency-based approach allowed to phenotype hundreds of populations per year, without the need for laborious strain isolation, thereby providing sufficient spatial coverage annually. The data were collected in a longterm network including diverse stakeholders, demonstrating the power of collective action for monitoring and managing fungicide resistance. Our analysis, based on field data, detected some original spatial cryptic structures and differences in the rate of increase of resistance frequencies between phenotypes and regions, making it possible to describe resistance status in France more precisely than the partial short-term analyses generally performed. This approach should improve regional resistance management. It offers a vision complementary to that provided by the theoretical models of resistance evolution generally proposed, and may also provide empirically validated parameter such as estimates of growth rates. Moreover, the three statistical models used here are not specific to Z. tritici and may be considered as generic tools, suitable for use in similar studies of other resistances. The value of this approach for describing resistance dynamics is now well established. Future studies should seek to explain the sources of variation in local resistance evolution, and to identify the extrinsic or intrinsic factors determining resistance dynamics, as these factors may be useful operational levers for resistance management. This could be achieved by developing the DYNAMIC MODEL further. Regional selective pressure (i.e. the regional use of fungicides) should be explored, together with demographics (e.g. the intensity of local epidemics) and local agronomic and environmental factors. Such improvement to the model would facilitate prediction of the structure of resistant populations year after year, while incorporating any context changes that could affect resistance evolution. Its outputs may help farmers to yearly adapt their spray programs while implementing sustainable local strategies before resistance reaches fixation in population. Indeed, according to the regional resistance dynamics and pesticide use, tailor-made strategies may be recommended, as mixture or alternation regimes may differ in efficacy locally. Like Hicks et al. 68 , in their long-term nationwide study of herbicide resistance, we conclude that there is a need to use "an evolutionarily informed approach in a proactive not reactive manner". [START_REF] Hardin | The tragedy of the commons[END_REF] with the exception of changes affecting the DYGYG motif. The L50S, D107V, N178S, S188N, S208T, S259F, N284H, H303Y, A311G, G312A, A410T, G412A, G476S, V490L, N504K, G510C and/or N513K variants are not related to azole resistance but are regularly reported in most genotypes, and would therefore greatly increase the total number of genotypes if reported. ¤ Nomenclature described in Huf et al. [START_REF] Grimmer | Fungicide resistance risk assessment based on traits associated with the rate of pathogen evolution[END_REF] . As more codons are considered in this nomenclature, several codes were possible to match with our genotype numbers.

TriS TriS G1 - D V Y A I Y G Y S 0 G2 A5 D V Y A I Y - Y S 1 TriLR TriR2 G3 A1 D V Y A I Y D Y S 1 TriLR TriR3 G4 A6 D V F A I Y G Y S 1 TriR4 G5 D1, D4,D5, E1, E14 D V Y A I - - Y S 1 G6 - D V Y A I H G Y S 1 G7 B8, D24 D V Y A I Y G S S 1 G8 A11, B9, D3 D V Y A I Y G H S 1 G9 E16, F3, F13 D C Y A I - - Y S 1 G10 B10 D C Y A I Y G H S 1 TriR5 G11 C7 D A Y A I Y G H S 2 G12 B2, C6 D A Y A I Y G S S 2 G13 C10, D8, E6, E11, E15, F5 D A Y A I - - Y S 2 TriMR TriR6 G14 B12, C1, C8, C16, D20, D21 D V Y A V Y G H S 2 X X G15 - D V Y A V Y G S S 2 G16 B4, C2 D V Y A V S G Y S 2 G17 B3, C4, D22 D V Y A V D G Y S 2 G18 - D V Y A V N G Y S 2 TriR7 G19 E17, E22, F2, F14 D V Y A V - - Y S 2 X X TriR8 G20 E20, F11, F16, G1, G4, H1, H2, H3, H5, H8 D V Y G V - - Y S 3 X TriHR TriR9 G21 C11, D7 D A Y A I Y G S T 3 - G22 - D C Y A I Y G H T 2 TriR10 G23 C12, D10 D A Y A V Y G H S 3 X X - G24 - D A Y A V Y G S S 3 X X - G25 E5 D A Y A V Y G H T 4 X X X TriR12 G26 E3 D A Y A V Y G S T 4 TriR11 G27 D18, E4 G A Y A V Y G H S 4 X X X - G28 F8, G5 G A Y A V Y G H T 5 X - G29 D13, E7, F4 D C Y A V Y G H T 3 - G30 - G V Y A V Y G H S 3 X X - G31 F7, G7 D A Y G V Y G S T 5 - G32 E9 G A Y A V S G Y S 4 - G33 F1 D V Y G V D G Y T 4 TriR9 G34 E23, F6, G3 D A Y A I - - Y T 3 TriRz G35 F12 G A Y A I - - Y S 3 - G36 - G C Y A I - - H S 3 - G37 G2, H9 D V Y A V - - Y T 3 - G38 - G V Y A V - - H S 4 X - G39 H4, I1 D A Y G V - - Y T 5 - G40 F9, H6, I3 D C Y G V - - Y T 4 - G41 - G C Y G V - - H T 6 X - G42 - G V Y G V - - H S 5
§ Description of additional resistance mechanisms possibly affecting strain susceptibility to azoles. Two types of insertion have been described in the promoter sequence of cyp51 and the shortest one has been shown to lead to cyp51 overexpression 6,7 . Three insertions (types 1-3) in the promoter sequence of mfs1 have been shown to induce overexpression, causing multidrug resistance (enhanced efflux from the membrane transporter MFS1, 8,9 ). 'X' indicates that the insertion was regularly found in the corresponding cyp51 background, as the three resistance mechanisms are independent. Some combinations may be missing. Estimates of interactions between region and year are not shown because the factor was not significant for all phenotypes. Estimates are given on a logit scale. Disregarding the value of the Year:Region interaction parameter, frequencies can be recovered (in %) by calculating the inverse logit of estimates, weighting by the estimated proportions of 0 and 100 (π0 and π100, respectively), and multiplying the value obtained by 100. For instance, the frequency estimated for StrR in 2011, in PCH at the T0 sampling time was: [0.02*0 + (1-0.02-0.29) *logit -1 (1.98+1.61-0.5-0.53)*)*(1-0.29-0.02)0.69 + 1*0.29*1 + 0*0.02]*] * 100 = 93%. Empty cells indicate that no data were available for estimation. Factors tests P-values are indicated in bold and italic. P-value thresholds: "•" (P < 0.1), "*" (P < 0.05), "**" (P < 0.025), "***" (P < 0.001). Regional adjustments of the initial national frequency were also estimated but are not shown here.

model{ # # Likelihood # for(i in 1:n){ y[i] ~ dbin(P[i], 100) # P[i] <-p0[i] * 1.0E-10 + pp[i] * p[i] + p100[i] * (1-1.0E-10) # p0[i] <-equals(inflated[i, phenotype[i]], 1) pp[i] <-equals(inflated[i, phenotype[i]], 2) p100[i] <-equals(inflated[i, phenotype[i]], 3) inflated[i, phenotype[i]] ~ dcat(pi[, phenotype[i]]) # p[i] <-exp(eta[i])/(1 + exp(eta[i])) # eta[i] ~ dnorm(mu[i], Phenotype_tau[phenotype[i]]) mu[i] <-Phenotype.YearXPhenotype.RegionXYearXPhenotype[region[i], year[i], phenotype[i]] + RegionXPhenotype[region[i], phenotype[i]] + Sampling.dateXPhenotype[sampling.date[i], phenotype[i]] } # # Inflated parameters # for(i in 1:n_phenotype){ pi[1:3, i] ~ ddirch(alpha[]) } # for(i in 1:3){ alpha[i] <-1/2 } # # Priors # for(i in 1:n_phenotype){ Phenotype[i] ~ dnorm(0, 1.0E-6) # for(j in 1:n_year){ Phenotype.YearXPhenotype[j, i] ~ dnorm(Phenotype[i], invVar.YearXPhenotype[i]) YearXPhenotype[j, i] <-Phenotype.YearXPhenotype[j, i] -Phenotype[i] # for(k in n_region){ Phenotype.YearXPhenotype.RegionXYearXPhenotype[k, j, i] ~ dnorm(Phenotype.YearXPhenotype[j, i], invVar.RegionXYearXPhenotype[k, i]) RegionXYearXPhenotype[k, j, i] <- Phenotype.YearXPhenotype.RegionXYearXPhenotype[k, j, i] -Phenotype.YearXPhenotype[j, i] } } # for(k in 1:n_region){ RegionXPhenotype[k, i] ~ dnorm(0, invVar.RegionXPhenotype[i]) } # Sampling.dateXPhenotype[1, i] ~ dnorm(0, 1.0E-6) # # # invVar.YearXPhenotype[i] <-1/(Sigma.YearXPhenotype[i]*Sigma.YearXPhenotype[i]) Sigma.YearXPhenotype[i] ~ dlnorm(0, 1) # invVar.RegionXPhenotype[i] <- 1/(Sigma.RegionXPhenotype[i]*Sigma.RegionXPhenotype[i]) Sigma.RegionXPhenotype[i] ~ dlnorm(0, 1) # for(k in 1:n_region){ invVar.RegionXYearXPhenotype[k, i] <-1/(Sigma.RegionXYearXPhenotype[k, i]*Sigma.RegionXYearXPhenotype[k, i]) Sigma.RegionXYearXPhenotype[k, i] ~ dlnorm(0, 1) } # Phenotype_tau[i] <-1/(Sigma.Phenotype_tau[i]*Sigma.Phenotype_tau[i]) Sigma.Phenotype_tau[i] ~ dlnorm(0, 1) } # # Constraints # for(i in 1:n_phenotype){ Sampling.dateXPhenotype[2, i] <-0 } } R script (JAGS): model{ # # Likelihood # for(i in 1:n){ y[i] ~ dbin(P[i], 100) # P[i] <-p0[i] * 1.0E-10 + pp[i] * p[i] + p100[i] * (1-1.0E-10) # p0[i] <-equals(inflated[i, phenotype[i]], 1) pp[i] <-equals(inflated[i, phenotype[i]], 2) p100[i] <-equals(inflated[i, phenotype[i]], 3) inflated[i, phenotype[i]] ~ dcat(pi[, phenotype[i]]) # p[i] <-exp(eta[i])/(1 + exp(eta[i])) # eta[i] ~ dnorm(mu[i], Phenotype_tau[phenotype[i]]) mu[i] <-Phenotype.RegionXPhenotype_T0[region[i],
Estimates are given on a logit scale. Initial frequency refers to frequencies observed in the column "First year" of Table 1. To recover frequencies (in %) the inverse logit of estimates must be calculated, weighted by the estimated proportion of 0 and 100 (π0 and π100, respectively), and the value obtained should then be multiplied by 100. For instance, for the estimated frequency of StrR in 2011 in PCH at the T0 sampling: [logit -1 (-1.04+(0.68+0.34)*7-0.58)*0.69 + 1*0.3 + 0*0.02]*100 = 98%. Empty cells indicate that no data were available for the estimation. P-value thresholds: "•" (P < 0.1), "*" (P < 0.05), "**" (P < 0.025), "***" (P < 0.001).

Figure 1 .

 1 Figure 1. Distribution of Performance network trials between 2004 and 2017 in mainland France. Dots represent trials and lines the regional structure before the reorganisation of French regions in 2015. Acronyms denote regions: ALS, Alsace; AUV, Auvergne; AQU, Aquitaine; BNO, Basse-Normandie; BOU, Bourgogne; BRE, Bretagne; CEN, Centre, CHA, Champagne-Ardennes; FCO, Franche-Comté; HNO,

Figure 2 .

 2 Figure 2. Changes in fungicide use and resistance frequency in Z. tritici populations in France.

Figure 3 .

 3 Figure 3. Apparent relative fitness estimates (increases in frequency of resistance) for different resistant phenotypes.

Figure 4 .

 4 Figure 4. Spatial partitioning over time for the StrR, TriR6 and TriR7-TriR8 phenotypes in France. Background maps with colour gradient are obtained by kriging methods, from real observations (points). The colours within the observation points indicate the true frequency observed in trials. The lines are obtained by the SPODT method (spatial oblique decision tree 53 ) and split space into two homogeneous blocks, one with lower frequencies than the other (or no resistance at all for the StrR phenotype). Left: background map and points representing the status of StrR resistance in 2006; colonization fronts from 2004 to 2007 are shown; for 2008, we show only the farthest trial from the 2007 front, because no more trials were found with frequencies close to 0% in this year. Middle: background map and points representing the status of TriR6 resistance in 2013; delimitation lines from 2006 to 2013 are shown. Right: background map and points representing the status of TriR7-TriR8 resistance in 2013; delimitation lines from 2006 to 2013 are shown.

Figure 5 .

 5 Figure 5. Distribution of mutations (insertion into the promoter of the MFS1 transporter gene) among MDR and non-MDR phenotypes in isolated Z. tritici strains (2008-2017; n=303).

Table 1 .

 1 The internal panels compiled, for each compound used for STB control, the amounts used in the region (in kg) and the corresponding deployed surfaces (in hectares). For a given fungicide, deployed hectares correspond to the number of sprays multiplied by the treated area, without regard to the dosage. For instance, 200 deployed hectares can mean either two applications of a single active ingredient over 100 ha, or one application of two different fungicides from the same mode of action at once over the same 100 ha. Panel data were available for all regions of France, from 1990 onwards. As an indicator of fungicide use, we calculated the national deployed hectares per mode of action (sum over all regions of the deployed surfaces for all the active ingredients sharing the same mode of action, used for STB control) divided by the area under wheat in France (from the AGRESTE Data used in the statistical analysis, extracted from the Performance database.

	Information about fungicide use collected via
	internal panels was provided by Bayer
	CropScience. Fungicide use was estimated
	from field surveys carried out on
	representative farms.

database: agreste.agriculture.gouv.fr). This indicator, denoted 𝐻𝐴𝐷 𝐶 (deployed hectares divided by cultivated hectares), represents the mean number of sprays of a mode of action in France for a given year. As the use of SDHIs is recommended to one spray per year in France, 𝐻𝐴𝐷 𝐶 for SDHIs represents the proportion of surfaces sprayed with this class of fungicides.

Table 2 .

 2 Frequencies and frequencies differences estimated from the ANOVA model for year, region and sampling date effects.

	Phenotype

Resistance ¤ Biochemical mode of action § Chemical structures § Active ingredient † Authorisation in France ‡ Codes Resistant phenotype abbreviation Target site resistance Non target site resistance Main processes Target site and mode of action Abbreviation Main classes Secondary classes Registration

  Legal authorisation of the active ingredient in France. Year the molecule was registered in France, year it was first used in the field and year it was removed from French registration (http://www.ephy.anses.fr).

	Imidazoles	Prochloraz	1980 1980?	-		[LR]
		Bromuconazole 1994 1995	-		Cyp51: one change (e.g.
		Cyproconazole	1987 1987	-		Y137F, Y459C/D/H,
		Difenoconazole 1988 2011	-		G460D/S, Y461H/S,
	First year of use 2016 2017 2006 2007 Benzovindiflupyr 2016 2017 Fluopyram Boscalid Bixafen 2011 2011 Penthiopyrad 2014 2014 Fluxapyroxad 2011 2011 Azoxystrobin 1996 1997 Picoxystrobin 1999 2000 Methoxy-carbamates Pyraclostrobin Benzamides Nicotinamides (syn. Pyridine-carboxamides) Pyrazole-carboxamides Methoxy-acrylates 2000 2000 Oximino-acetates Kresoxim-methyl 1996 1997 2013 Removal ---------Trifloxystrobin 2001 2002 -Oximino-acetamides Dimoxystrobin 2009 2010 -Fluoxastrobin 2006 2006 -Epoxiconazole 1992 1992 -Triazoles Fenbuconazole 1991 1992 -Fluquinconazole 1997 1998 2016 Flusilazole 1985 1986 2013 Flutriafol 1983 1984 2015 Hexaconazole 1990 1990 2007 Metconazole 1993 1994 -Propiconazole 1980 1980? -Tebuconazole 1988 1989 -Tetraconazole 1991 1991? -Triadimenol 1987 1988 2014 Triazolinethiones Prothioconazole 2006 2006 -Benzimidazoles Carbendazim 1972 1972? 2009 Thiophanates Thiophanate-methyl 1973 1980? -Akylene-bis-dithiocarbamates Mancozeb 1980 1980? -Maneb 1988 1988? -Chloronitriles Phthalonitriles SDHI Carboxamides QoI (or QoI-P Synthetic strobilurins and analogs DMI or SBI-DM Heterocyclic compounds MBC Benzimidazole s and precursors MSI Dithio-carbamates Chlorothalonil 1980 1980? -Halogenated alkyl-thio coumpounds Phthalamides Folpet 2003 2003 2009 U-W10 M4 R4P § FRAC † U-A2a C2/7 U-A5a C3/11 U-E2 G1/3 U-K2b B1/1 U-W11 M3 U-W14 M5 § According to the R4P classification (http://www.r4p-inra.fr). Mitochondrial respiration Complex II or succinate dehydrogenase. Binding site of ubiquinone involving the subunits SdhB, SdhC and SdhD. Complex III or cytochrome bc1. Binding site of ubiquinone to cytochrome b at the "o" center (proton output) in the heme bl proximal domain Sterol biosynthesis Binding to sterol 14α-demethylase (syn. CYP51) Cell division: microtubules of the mitotic spindle Binding to β-tubulin Multis-site activity Fungicide acting non specifically on multiple targets, especially respiratory enzymes ? Approximate date. NC: Not concerned. † According to the FRAC classification (http://www.frac.info). Main phenotype class † Historical phenotype † Genotype number cyp51 codon ‡ Genotype in nomenclature ¤ 134 136 137 379 381 459 460 461 524	CarR StrR TriR BenR NC Number of changes	+ + + + NC	Target alteration [LR] SdhB: N225I, R265P, T268A/I SdhC: T79N, W80S, N86A/S R151M/S/T, I261S, V166M, T168R SdhD: R47W, I50F, M114F [MR-HR] SdhC: T79I, H152R [LR] Cytb: F129L [HR] Cytb: G143A ΔY459/G460) or two changes (V136A+Y461H/S or ΔY459/G460) [MR] Cyp51: two changes including I381V+Y459S/D/N, Y461H/S or ΔY459/G460) [HR] Cyp51: three to six changes including those reported for TriLR and TriMR phenotypes, combined with D107V, D134G and/or S524T. [HR] Tub2: E198A NC Insert detection § Target overexpression --+ [LR when alone] Cyp51 promoter: 120 bp or 1000 bp inserts -NC cyp51 cyp51 mfs1 mfs1 expression Compensation -+ [LR] AOX over---NC mfs1 promoter: Enhanced efflux (MDR) + [LR when alone] MFS1 types I, II or III inserts + [LR when alone] MFS1 promoter: types I, II or III inserts + [LR when alone] MFS1 promoter: types I, II or III inserts -NC 120 bp 1000 bp type 1 type 2 type 3

‡

  X † Nomenclature described in Leroux et al. 1 and Leroux & Walker 2 .‡ Numbers indicate codons affected by mutation or deletion within the cyp51 sequences of the strains of the collection. Letters indicate amino acids according to standard nomenclature. For the sake of simplicity, only codons relevant for azole resistance are reported, according to Cools et al.3 and Cools et al.

  Leroux P, Albertini C, Gautier A, Gredt M, and Walker A-S, Mutations in the cyp51 gene correlated with changes in sensitivity to sterol 14α-demethylation inhibitors in field isolates of

		Mycosphaerella graminicola, Pest Manag Sci 63:688-698 (2007).
	2	Leroux P and Walker A-S, Multiple mechanisms account for resistance to sterol 14α-
		demethylation inhibitors in field isolates of Mycosphaerella graminicola, Pest Manag Sci
		67:44-59 (2011).
	3	Cools H, Parker J, Kelly D, Lucas J, Fraaije B, and Kelly S, Heterologous expression of mutated
		eburicol 14α-demethylase (CYP51) proteins of Mycosphaerella graminicola to assess effects
		on azole fungicide sensitivity and intrinsic protein function, Appl Environ Microbiol 76:2866-
		2872 (2010).
	4	Cools HJ, Mullins JG, Fraaije BA, Parker JE, Kelly DE, Lucas JA, et al., Impact of recently emerged
		sterol 14α-demethylase (CYP51) variants of Mycosphaerella graminicola on azole fungicide
		sensitivity, Appl Environ Microbiol 77:3830-3837 (2011).
	5	Huf A, Rahfus A, Lorenz K, Bryson R, Voegele R, and Stammler G, Proposal for a new
		nomenclature for CYP 51 haplotypes in Zymoseptoria tritici and analysis of their distribution in
		Europe, Plant Pathol (2018).
	6	Chassot C, Hugelshofer U, Sierotzki H, Gisi U, and others, Sensitivity of CYP51 genotypes to DMI
		fungicides in Mycosphaerella graminicola, Modern fungicides and antifungal compounds V:
		15th International Reinhardsbrunn Symposium, Friedrichroda, Germany, May 6-10, 2007.
		Deutsche Phytomedizinische Gesellschaft eV Verlag, 129-136, (2008).

7 Cools HJ, Bayon C, Atkins S, Lucas JA, and Fraaije BA, Overexpression of the sterol 14αdemethylase gene (MgCYP51) in Mycosphaerella graminicola isolates confers a novel azole fungicide sensitivity phenotype, Pest Manag Sci 68:1034-1040 (2012). 8 Omrane S, Audeon C, Ignace A, Duplaix C, Sghyer H, Aouini L, et al., Multi-drug-resistance (MDR) in Septoria leaf blotch, Resist 2015 67, 44-59 (2015). 9 Omrane S, Audéon C, Ignace A, Duplaix C, Aouini L, Kema G, et al., Plasticity of the MFS1 promoter leads to multidrug resistance in the wheat pathogen Zymoseptoria tritici, mSphere 2:e00393-17 (2017). Supplementary Information 3.1. ANOVA MODEL Model: R Script (JAGS):

  Estimates for the ANOVA MODEL, for year, region and sampling date effects.

	Supplementary Information 4.
	phenotype[i]] +
	Phenotype.RegionXPhenotype[region[i], phenotype[i]]*time[i] +
	Sampling.dateXPhenotype[sampling.date[i], phenotype[i]]
	}
	#
	# Inflated parameters
	#
	for(i in 1:n_phenotype){
	pi[1:3, i] ~ ddirch(alpha[])
	}
	#
	for(i in 1:3){
	alpha[i] <-1/2
	}
	#
	# Priors
	#
	for(i in 1:n_phenotype){
	Phenotype_T0[i] ~ dnorm(0, 1.0E-6)
	#
	Phenotype[i] ~ dnorm(0, 1.0E-6)
	#
	for(j in 1:n_region){
	Phenotype.RegionXPhenotype_T0[j, i] ~ dnorm(Phenotype_T0[i],
	invVar.RegionXPhenotype_T0[i])
	RegionXPhenotype_T0[j, i] <-Phenotype.RegionXPhenotype_T0[j, i] -
	Phenotype_T0[i]
	#
	Phenotype.RegionXPhenotype[j, i] ~ dnorm(Phenotype[i],
	invVar.RegionXPhenotype[i])
	RegionXPhenotype[j, i] <-Phenotype.RegionXPhenotype[j, i] -
	Phenotype[i]
	}
	#
	Sampling.dateXPhenotype[1, i] ~ dnorm(0, 1.0E-6)
	#
	#
	#
	invVar.RegionXPhenotype_T0[i] <-
	1/(Sigma.RegionXPhenotype_T0[i]*Sigma.RegionXPhenotype_T0[i])
	Sigma.RegionXPhenotype_T0[i] ~ dlnorm(0, 1)
	#
	invVar.RegionXPhenotype[i] <-
	1/(Sigma.RegionXPhenotype[i]*Sigma.RegionXPhenotype[i])
	Sigma.RegionXPhenotype[i] ~ dlnorm(0, 1)
	#
	Phenotype_tau[i] <-1/(Sigma.Phenotype_tau[i]*Sigma.Phenotype_tau[i])
	Sigma.Phenotype_tau[i] ~ dlnorm(0, 1)
	}
	#
	# Constraints
	#
	for(i in 1:n_phenotype){
	Sampling.dateXPhenotype[2, i] <-0
	}
	}

  Supplementary Information 5.Estimates for the DYNAMIC MODEL.

							Phenotype
	Factor	Parameter	BenR	StrR	TriLR	TriR2-TriR4	TriR5	TriMR	TriR6	TriR7-TriR8	TriHR	MDR
	Inter-annual national mean	1.78***	1.98* -1.97*** -2.1*** -2.52*** 1.19***	0.31	-1.32*** -1.87*** -1.99***
	Year		0.894	0***	0***	0.065 •	0***	0***	0***	0***	0.034*	0.972
	2004		-3.04**
	2005		-1.93*	1.7***			-0.83**
	2006		-0.89	0.99*	0.31*	0.45	-0.25	-0.15	-0.07
	2007	-0.2	-0.57	0.76 •	0.27 •	0.72*	-0.05	-0.15	0.16
	2008	-0.14	-0.56	0.76 •	0.31*	0.55*	-0.03	-0.15	0.09
	2009	0.14	-0.34	0.28	-0.03	0.47 •	0.42*	0.5*	-0.35
	2010	0.03	0.94	0.13	-0.2	0.23	0.4*	0.5*	-0.36	-1.43**
	2011	0	1.61 •	-0.06	-0.17	-0.04	0.5*	0.56*	-0.19	-0.58
	2012	0.24	4.84***	-0.04	-0.04	-0.29	0.17	0.35	-0.74**	-0.42
	2013	0.16		0.18	0.18	-0.1	0.04	-0.55*	0.52*	0.05	-0.18
	2014	0.03		-0.54	-0.32	-0.71 •	0.15	0.45	-0.18	0.38	-0.08
	2015	-0.18		-1.46**	-0.3	-1.26**	0.37	0.39	0.34	0.13	0.22
	2016	0.03		-2.87***			-0.55*	-0.14	-0.54	1.07*	0.1
	2017	-0.09					-0.41 -1.57*** 1.37***	0.8 •	-0.05
	Region		0***	0.224 0.004**	1	0.018*	0.99	0.01*	0***	1	1
	NPC	0.58*	1.06*	-0.49*	-0.06	-0.41	0.17	0.54*	-0.64*	0.34	-0.02
	PIC		0.56**	0.47	-0.59*	-0.18	-0.6*	0.16	0.4 •	-0.52 •	0.27	0
	HNO	0.61 •	0.79*	-0.3	-0.06	0.02	0.15	0.44 •	-0.43	0.01
	BNO	0.18	0.23	0.08	0.09	-0.07	-0.14	0.18	-0.66*	0.39
	LOR	0.35	-0.02	0.15	0	-0.17	-0.03	0.34	-0.46	-0.19
	CHA	0.25	0.19	-0.05	-0.04	-0.19	0.05	0.29	-0.45 •	0.15	-0.03
	IDF		0.27	0.43	-0.1	-0.03	0.03	0.07	0.53*	-0.46	-0.27	-0.05
	BRE	-0.67**	0.21	-0.27	-0.02	-0.33 •	0.07	-0.23	0.42	0.12
	CEN	0.03	0.4	0.34*	0.02	0.38 •	-0.07	-0.03	-0.16	-0.35	0.01
	PDL	-0.48 •	-1.08*	0.28	0.02	0.6*	-0.17	-0.37 •	0.31	-0.05	0.09
	BOU	0.05	-0.22	0.17	0.05	0.15	-0.02	-0.04	-0.14
	PCH	-0.2	-0.5	0.57**	-0.04	0.66**	-0.22	-0.45	0.63*	-0.13
	AUV	0.04	-0.09	0.23	0.09	0.17	-0.17	-0.08	-0.22	0.15
	RAL	-0.32	-0.1	0.14	0.04	-0.03	-0.03	0.18	-0.07	-0.33
	AQU	-0.91** -1.19 •	0.06	0.07	0.1	0.03	-0.64 •	0.94*	-0.12
	MPY	-0.33	-0.53	-0.21	0.06	-0.27	0.18	-1.1 •	1.88**
	Year:Region		1	0.661	1	1	1	1	1	1	1	1
	Sampling date TNT-T0	0.04	0.53*** -0.21** -0.25**	-0.06	0.2**	0.39***	-0.11	0.61
	π0		0	0.02	0.1	0.3	0.4	0	0.05	0.17	0.47	0.61
	π100	0.3	0.29	0	0	0	0.07	0.03	0	0	0
	R²		0.45	0.84	0.59	0.41	0.53	0.42	0.64	0.58	0.47	0.6

ACKNOWLEDGEMENTS

This work was funded by the INRA SMaCH metaprogram (FONDU project) and by Bayer CropScience (PhD scholarship of M. Garnault; CIFRE no. 2016-0695). The authors warmly thank all partners of the Performance network, for providing samples since 2004, and Arvalis-Institut du Végétal for supervising field work. The authors thank the INRA MIGALE bioinformatics platform (http://migale.jouy.inra.fr) for providing computational resources, Laurent Thibault (Anses) for proving information about fungicide registration, and Robert Canalès, Constance Tuffet and Clarisse Payet (Bayer CropScience) for providing data on fungicide use in France and for helpful discussions.

Fungicides

Supplementary Information 1. Fungicides used to control STB in France and associated resistance phenomena. Description and resistance mechanisms are detailed in the following references: BenR 1,2 ; TriR [START_REF] Barres | Trends and challenges in pesticide resistance detection[END_REF][START_REF] Hardin | The tragedy of the commons[END_REF] ; StrR [START_REF] Delp | Fungicide resistance: definitions and use of terms[END_REF][START_REF] Grimmer | Fungicide resistance risk assessment based on traits associated with the rate of pathogen evolution[END_REF][6][7] ; CarR 8-10 ; MDR [START_REF] Barres | Trends and challenges in pesticide resistance detection[END_REF][11][12][13] . Amino-acid changes were described in the corresponding publications. Supplementary Information 2. cyp51 genotypes in a collection of pure strains. Strains were isolated by INRA between 1969 and 2017 (n=314). This list is not exhaustive. Other genotypes may exist in populations and have been reported in other studies.

Supplementary Information 3.2. DYNAMIC MODEL