Patterns of thermal adaptation in a globally distributed plant pathogen: Local diversity and plasticity reveal two‐tier dynamics

Abstract Plant pathogen populations inhabit patchy environments with contrasting, variable thermal conditions. We investigated the diversity of thermal responses in populations sampled over contrasting spatiotemporal scales, to improve our understanding of their dynamics of adaptation to local conditions. Samples of natural populations of the wheat pathogen Zymoseptoria tritici were collected from sites within the Euro‐Mediterranean region subject to a broad range of climatic conditions. We tested for local adaptation, by accounting for the diversity of responses at the individual and population levels on the basis of key thermal performance curve parameters and “thermotype” (groups of individuals with similar thermal responses) composition. The characterization of phenotypic responses and genotypic structure revealed the following: (i) a high degree of individual plasticity and variation in sensitivity to temperature conditions across spatiotemporal scales and populations; and (ii) geographic variation in thermal response among populations, with major alterations due to seasonal patterns over the wheat‐growing season. The seasonal shifts in functional composition suggest that populations are locally structured by selection, contributing to adaptation patterns. Further studies combining selection experiments and modeling are required to determine how functional group selection drives population dynamics and adaptive potential in response to thermal heterogeneity.


| INTRODUC TI ON
Environmental heterogeneity (Li & Reynolds, 1995) is regarded as one of the most important elements driving the emergence and maintenance of genetic variation within populations (Hedrick, 1986;Hughes et al., 2008;Levins, 1974;Ravigné et al., 2009) as it dictates physiological responses (Cavieres & Sabat, 2008) and can drive the emergence of local adaptation patterns (Nuismer & Gandon, 2008;Thompson, 2005). Gathering information about the way a given community, species, or population copes with this environmental heterogeneity is crucial for the understanding and prediction of its distribution and responses to current and future environmental changes (Austin, 2007).
The adequate capture of eco-evolutionary responses requires an integration of physiological variation across biological (individual, group, population, species) and spatiotemporal (seasonal, geographic) scales, given the significant implications of this variation for dynamics (Saloniemi, 1993;Schreiber et al., 2011;Vindenes et al., 2008). It is therefore important to go beyond summarizing diversity through average trait values Violle et al., 2012), and to account for the individual specialization of phenotypic responses by taking into account both phenotypic plasticity (withinindividual differences; Pigliucci, 2001) and interindividual variation (between-individual differences; Dall et al., 2012).
The ecological concept of "reaction norm," describing the set of phenotypes generated by a given genotype in different environments (Schlichting & Pigliucci, 1998), is particularly effective as a tool for accounting for individual specialization (Araújo et al., 2008;Bolnick et al., 2002). Most of the comparisons of the thermal sensitivity of a given phenotypic trait across individuals under different environmental conditions have been conducted to date on reaction norm descriptors (e.g., comparisons of mean phenotypic differences and cardinal temperatures; Gibert et al., 1998) or degree of plasticity (e.g., regression for linear reaction norms or horizontal [warmercolder], vertical [faster-slower], and shape [generalist-specialist] shifts of non-linear reaction norms; Izem & Kingsolver, 2005;Martin et al., 2011;van de Pol, 2012). Such approaches have proved highly valuable, but may not be suitable for decomposing the overall variation or distinguishing differential responses among populations (Bulté & Blouin-Demers, 2006) or including intra-and interindividual sources of error (Angilletta, 2006) in ANOVA and random regression approaches (Gilchrist, 1995;Lynch & Gabriel, 1987).
One possible complementary approach to the description of variation between reaction norms involves the use of functional ecology to describe significant variations in the degree of individual specialization within populations and species (Garnier & Navas, 2012). The idea is to translate reaction norms by grouping individual reaction norms into "functional groups" (Gitay & Noble, 1997;Violle et al., 2007). Each of these functional groups responds to the environment in its own way (e.g., low-or high-performance specialists), according to a classification system that is not predetermined (i.e., constrained modes of variation). This approach accounts more effectively for patterns of variation in phenotypic plasticity, through the characterization of three functional components: richness, evenness, and divergence (Mason et al., 2005).
This approach is particularly useful for deciphering variation in continuous reaction norms describing performance as a function of temperature (thermal performance curves or TPC; Huey & Stevenson, 1979), and for documenting patterns of thermal adaptation to prevailing local conditions (Kawecki & Ebert, 2004) across a range of environments (e.g., Mitchell & Lampert, 2000). These patterns play an important role in the case of microorganisms impacting ecosystems, human health, and food security (Fisher et al., 2012) as local adaptation to temperature conditions governs their geographic distribution, phenology, and abundance (Kraemer & Boynton, 2017). This results in impacting the expansion ranges of plant pathogens (e.g., Milus et al., 2009;Robin et al., 2017), as well as the onset and severity of disease epidemics (e.g., Ferrandino, 2012).
In plant pathogens, summarizing the individual variance of aggressiveness traits as population-scale averages (problematic use of single mean species values; Suffert & Thompson, 2018) or phenotyping individuals under a limited set of temperatures when considering variances (generally about three temperatures in thermal biology studies; Dell et al., 2013;Low-Décarie et al., 2017) has provided useful information about species distribution. This has made it possible to detect signatures of interindividual variation and adaptation within species and populations (Milus et al., 2006). However, this information cannot be used to infer selection driving population dynamics (Lavergne et al., 2010) or to assess the relevant scales of functional diversity (Martiny et al., 2011;Woodcock et al., 2006). This study explored the extent of variation in thermal responses of a globally distributed wheat pathogen across space (geographic range) and time (local seasonal dynamics), and uncovered the role that adaptation to local environmental conditions (dynamic evolutionary process) plays in generating this diversity. The analysis of the plasticity and variation of thermal sensitivity across individuals, populations, and scales was conducted in the case of Zymoseptoria tritici (formerly Mycosphaerella graminicola; Steinberg, 2015), the causal agent of one of the most economically important wheat diseases (Septoria tritici blotch or STB; Dean et al., 2012;Fones & Gurr, 2015). Besides its agronomic relevance, we chose to study this fungal pathogen as its aggressiveness traits are empirically known to be temperature-sensitive (Lovell et al., 2004;Shaw, 1990) and to display interindividual variation (Bernard et al., 2013;Boixel et al., 2019). The duality of the reproduction modes-asexual and sexual, which both contribute to the local level of genetic structure (Singh et al., 2021;Suffert & Sache, 2011)-makes this epidemiological model particularly interesting: (i) Sexual lineages maintain and increase genetic diversity in pathogen populations, through sexual spores that are wind-dispersed over long distances from wheat residues at the end of each growing season; and (ii) clonal lineages (asexual reproduction) occur within a single field during the course of an epidemic, through asexual spores rain-dispersed over short distances. Furthermore, Z. tritici populations present signatures of adaptation to a wide range of contrasted environments over space (globally distributed pathogen across wheat-growing areas worldwide; Zhan & McDonald, 2011) and time (covering seasonal changes, e.g., from late autumn to early summer in Europe; Suffert et al., 2015). Drawing on previous local adaptation studies conducted by Zhan and McDonald (2011) and Suffert et al. (2015), we designed a sampling scheme to grasp the levels of functional diversity shaping responses of Z. tritici populations to contrasted environments.

| Sampling survey design
Samples were collected from 12 Z. tritici populations for the exploration of spatial and temporal components of thermal adaptation (one population being a sample of the complete group of individuals occupying a given wheat plot at a spatiotemporal location) (Figure 1-Step 1). Spatial variation was investigated for 8 populations sampled within the Euro-Mediterranean region (see detailed sampling information of the geographic scale in Table 1) representative of the contrasting climatic conditions over this large geographic area (covering three-Cfb, Csa, and Dfb-out of seven Köppen-Geiger climate zones in which Z. tritici is reported as a notable pathogen; Figure S1). One of these sites (Grignon, France) was selected for a comparison of the thermal responses of two pairs of winter and spring subpopulations sampled from neighboring fields, to capture seasonal dynamics over a wheat-growing season (i.e., over the course of an annual epidemic; see the 4 populations of the seasonal scale in Table 1 and Figure S2). These pairs of subpopulations were subject to seasonal variation from November to February and from March to June, respectively. For each of the 12 populations, we collected 25-30 isolates at random from wheat leaves with STB symptoms, which were placed on a wet filter paper in moist chambers to promote the extrusion of Z. tritici cirrhi. On each leaf, one cirrhus from a single pycnidium was retrieved for isolation in pure culture (see Methods S1 for more details). After two subculturing for obtaining pure single-spore strains, Z. tritici spore suspensions were stored at −80°C in cryotubes, in a 1:1 glycerolwater mixture. Prior to thermal phenotyping, strains were subcultured once just after their thawing. The time elapsed between strain isolation and phenotyping experiments ranged from 1 (Euro-Mediterranean geographic populations) to 5 (French seasonal subpopulations) years. These conditions reduced the potential effects of previous environmental acclimation (for instance, via transgenerational plasticity or epigenetics), although such effects cannot be completely excluded. The strains were later confirmed to be genetically unique strains based on neutral genetic markers. We chose to consider 25 or 30 strains (i.e., individuals) per population instead of the minimum level of 15 identified on the basis of a rarefaction analysis ( Figure S3) for estimating the diversity of thermal responses in Z. tritici with sufficient power, accuracy, and precision (Dale & Fortin, 2014).

| Phenotypic variations in thermal responses
Thermal responses were phenotyped by determining the in vitro growth rates of the strains in liquid glucose peptone medium (14.3 g L −1 dextrose, 7.1 g L −1 bactopeptone, and 1.4 g L −1 yeast extract) over a 4-day period at 12 constant temperatures ranging from 6.5 to 33. 5°C (6.5, 9.5, 11.5, 14.5, 17.5, 20.0, 22.5, 24.5, 26.5, 28.5, 30.5, and 33.5°C;Boixel et al., 2019) (Figure 1-Step 2). The growth rate μ of each strain at each temperature (n = 8; independent replicates) was calculated according to the standardized specific experimental framework developed by Boixel et al. (2019), which has been validated to be representative of in planta responses with respect to discrimination between "cold-and warm-adapted" individuals. Thermal performance curves (TPCs) describing in vitro growth rate as a function of temperature were established by fitting a quadratic function to the temperature-growth rate (or performance P) estimates for each strain: P(T) = P max + Curv(T − T opt ) 2 where Curv is a shape parameter (see Table S1 for more information on the selection process of the model leading to the highest accuracy of performance estimates over the mid-temperature range). The key properties of TPCs were estimated through thermal traits commonly used to compare thermal sensitivities (Angilletta, 2006;Kingsolver, 2004). We have retained three parameters to describe the shape of these TPCs and quantify their characteristics: first, maximum performance (P max ), which informs on TPC height F I G U R E 1 Overview of the methodology for characterizing diversity and adaptive patterns of thermal responses in the sampled Zymoseptoria tritici populations. (Step 1) Twelve populations, each composed of either 25 or 30 strains, were collected from diseased leaves in different spatiotemporal locations (8 Euro-Mediterranean populations collected along geographic thermal gradients and 4 French seasonal subpopulations) with the corresponding mesoclimatic conditions (temperature data). All strains were (Step 2) phenotyped in an in vitro growth experiment conducted over a range of 12 temperatures to capture thermal performance curves from growth kinetics involving  All individuals from a given location were collected from a single pure stand field of a bread wheat cultivar susceptible to STB, except for the TN population, which was sampled from a durum wheat cultivar. f Z. tritici populations were sampled during the 2015-2016 wheat-growing season, except for the IS population, which was sampled during the 2016-2017 growing season.
("vertical shift" modes of variation); second, thermal optimum (T opt ), which informs on TPC position at the peak performance ("horizontal shift" modes of variation); and third, thermal performance breadth (temperature range over which performance exceeds 80% of P max ; TPB 80 ), which informs on the sensitivity of the response to temperature change around T opt ("width shift" modes of variation).
The estimates of the minimum and maximum temperatures (T min and T max ), which define limits of growth, were not retained for further analysis as they fell outside the range of temperatures tested.
Differences in thermal responses were assessed in two successive ways: (i) differences in the range and mean values of P max , T opt , and TPB 80 , assessed with parametric or nonparametric (depending on whether the assumptions of normality and homoskedasticity were verified) statistical tests for comparing variances and means; and (ii) typological comparisons grouping together TPCs with similar thermal characteristics (functional thermal groups, referred to hereafter as "thermotypes") based on a K-means clustering procedure applied to the covariation of P max , T opt , and TPB 80 for all TPCs (Methods S2). This further analysis of TPCs in thermotypes allowed to categorize individuals into five classes of horizontal position of the curve: "highly cold-adapted" (CA + ), "cold-adapted" (CA), "intermediate" (-), "warm-adapted" (WA), and "highly warm-adapted" (WA + ) that perform better at lower, low, median, high, and higher temperatures, respectively. When in quotation marks here and hereafter, the term "adapted" refers to this higher performance at specific temperature ranges (e.g., cold or warm environment performers) and, at this point, not directly to an assumption about the adaptation to the local environment in which they have been collected. A comparison of the distribution patterns in thermotypes across populations and scales was conducted to detect phenotypic differentiation based on chi-squared tests on the observed frequency distribution of thermotypes.

| Neutral genetic variation and population differentiation
To assess population genetic differentiation, the 350 individuals composing the 12 Z. tritici populations were genotyped for 12 neutral genetic markers on DNA extracted from 50 mg of fresh fungal material from 5-day cultures, following SSR amplification and sequencing in one multiplex PCR sample, and allele size annotation (Gautier et al., 2014; Methods S3a) ( Figure 1-Step 3).
Population structure was inferred with a Bayesian clustering approach under an admixture and correlated allele frequency model implemented in STRUCTURE (Pritchard et al., 2000). The degree and significance of genetic variability within a population (genetic diversity and allele richness) and differentiation between populations (pairwise estimates of Weir and Cockerham's F-statistic-F ST -and hierarchical analyses of molecular variance-AMOVA) were evaluated with random allelic permutation procedures in GENETIX (Belkhir, 2004) and Arlequin (Excoffier & Lischer, 2010) software (Methods S3b-d).

| Characterization of local climates
Air temperature data for the closest weather stations within a mean  Figure S4).

| Testing for signatures of local adaptation
Two steps were taken to detect genetic and phenotypic signatures of local adaptation underlying the observed differentiation between populations ( Figure 1-Step 5). First, the degree of genetic differentiation for the set of neutral markers (F ST index; Weir & Cockerham, 1984) was compared with that for phenotypic traits (P ST index; Leinonen et al., 2006). This made it possible to infer departures from neutral expectations (Merilä & Crnokrak, 2001), to determine whether thermal traits were under selection rather than subject to genetic drift (Brommer, 2011).  Figure  S5). A five-level scale was defined to summarize the overall difference in low and high values of P max (low-vs. high-performance strains); T opt (cold-vs. warm-adapted strain); and TPB 80 (specialist vs. generalist strain): statistically significant (1) much lower, (2) lower, (3) no deviation, (4) higher, and (5) much higher value, relative to the overall mean of each parameter over the whole data set. The indicated T opt , TPB 80 , and P max values correspond to the "barycenter" of each thermotype. (b, c, d) Three common documented shifts in thermal biology studies were identified: (b) a horizontal shift with variations in the position of TPCs along the temperature axis distinguishing "cold-adapted" vs. "warm-adapted" thermotypes; a horizontal stretch distinguishing "generalist" vs. "specialist" thermotypes (c) without or (d) with tradeoffs between P max and TPB 80 (TPC axes: P: performance; T: temperature). (e) Scatter plot highlighting a trade-off between P max and TPB 80 . P max is generally negatively related to TPB 80 except for a group of TPCs with both high TPB 80 and P max (green points surrounded by a rectangle). The regression is displayed as a solid line, with its 95% confidence interval as a shaded area, together with Pearson's correlation coefficient R and its p-value p the three thermal traits at the seasonal scale (Figure 2b-

| A reading grid for functional diversity in individual thermal responses
TPCs were classified into thermotypes with similar thermal responses thermotypes (Th1 to Th13; Figure S5), for which relative degrees of temperature specialization were described in terms of the T opt ("coldvs. warm-adapted"), TPB 80 ("specialist vs. generalist"), and P max ("lowvs. high-performer") dimensions (Figure 3a). These thermotypes illustrated two commonly documented non-exclusive shifts in TPC along thermal gradients: a horizontal shift (low-temperature vs. hightemperature generalists or low-temperature vs. high-temperature specialists; e.g., Th1 vs. Th13 in Figure 3b) and a generalist-specialist shift without (Th8 vs. Th9 in Figure 3c) or with (Th1 vs. Th3 or Th11 vs. Th13 in Figure 3d) trade-offs between P max and TPB 80 (i.e., when one cannot increase without a decrease in the other). Indeed, regression analysis revealed a significant negative correlation between P max and TPB 80 across all individuals (Pearson's correlation coefficient: and Th10). The distinguishing features of these four thermotypes were their average behavior with respect to T opt (Th5, Th6, Th7), TPB 80 (Th5, Th10), and P max (Th7, Th10).

Mediterranean populations
For population-level TPCs, significant variation was observed for thermal trait means for T opt (Kruskal-Wallis test, p < .01) and TPB 80 (Kruskal-Wallis test, p < .01), but not for P max (Kruskal-Wallis test, p = .09), for which no population differentiation was detected (Table 2). P max values may have been constrained by the upper detection thresholds for optical density (potential saturation of absorbance measurements for individuals with "extreme performance phenotypes"). There was a two-degree difference in T opt between the IS population and the 7 other populations ( Table 2) Table 2). These two populations had opposite patterns in terms of their respective proportions of thermal specialists and generalists (Figure 4b and Figure S6b). More broadly, the individuals with the greatest thermal breadth (G + , Th1-Th13) were less abundant in Cfb populations (DK-FR-IR), which were characterized by a higher proportion of more highly specialist individuals (S + , Th4, and Th9) than the average (accounting for 10% of the total chi-squared score; Figure S6d).

| Seasonal phenotypic shifts within local populations
Spring subpopulations (SPR1 and SPR2) had a higher thermal opti-

Conversely, WIN populations had a higher proportion of individuals
performing better under colder conditions (CA + ; Figure S7).

| Signatures of local adaptation to mean annual temperature conditions
Neutral molecular markers revealed that all strains were genetically different. We observed no difference in the genetic structure

F I G U R E 6
Functional diversity in thermal responses between the 12 Zymoseptoria tritici populations. Geographic (in bold) and seasonal (in standard text) populations are situated along: (i) a scale of increasing degree of adaptation to warm conditions (y-axis) discriminating colder-and warmer-adapted populations (logarithm of the ratio of the total of "warm-adapted" individuals-WA and WA + -to the total of "cold-adapted" individuals-CA and CA + ); (ii) a scale of thermal breadth continuum (x-axis) discriminating more specialist and more generalist populations (logarithm of the ratio of the total number of generalist individuals-G and G + -to the total number of specialist individuals-S and S + ) normals) indicated that the mean thermal optimum of geographic populations increased with mean annual temperature (Figure 7a).
The level of cold adaptation of these populations (measured as the ratio of highly "cold-adapted" to highly "warm-adapted" strains) was negatively and significantly correlated with the same climatic variable ( Figure 7b).

| Thermal phenotyping of Zymoseptoria tritici strains beyond the usual tests of "thermal sensitivity"
As several other studies on thermal phenotyping (Birgander et al., 2018;Paisley et al., 2005;Robin et al., 2017;Stefansson et al., 2013;Zhan & McDonald, 2011), the high-throughput method used in . As such, this framework enables to detect differences in thermal sensitivity between isolates (whatever the physiological bases that underpin these differences) and to go beyond the usual tests of "thermal sensitivity" based on two temperatures, which can be misleading due to the non-linearity of reaction norms (Angilletta, 2009). An advantage of the in vitro approach is that it enables largescale investigations while alleviating major confounding factors (e.g., cross-effect between host resistance and temperature adaptation; Pariaud et al., 2009). It should be mentioned, however, that, in contrast to in vitro responses, in planta processes exhibit narrower temperature responses and shifts to lower thermal optima (Boixel et al., 2019;Chaloner et al., 2020). Despite this in planta restriction of temperature niche breadth, a ranking consistency of thermal sensitivity between "cold-and warm-adapted" strains, consistent with the concept of phenotypic integration (Pigliucci, 2003), has been reported in previous studies (e.g., Paisley et al., 2005), notably in the case of Z. tritici (Boixel et al., 2019;Zhan et al., 2016). This alteration of thermal responses related to disease in planta may be due to suboptimal resource conditions (e.g., interaction with the host plant, stress responses, and nutrient restriction) compared with growth in axenic culture (Chaloner et al., 2020). Pearson's correlation coefficient: R = .98, p < .01). Population differentiation in T opt relative to neutral genetic differentiation is indicated by P ST and F ST values. (b) Relationship between cold adaptation level, defined as CA + /WA + (ratio of the number of "highly cold-adapted" to "highly warm-adapted" thermotypes), and the mean annual temperature at the sampling site. Linear dependence between these pairs of variables is indicated by the regression line (solid line and its 95% confidence interval, shown as a shaded area), Pearson's correlation coefficient R, and its associated p-value p (see Figure S1 for a description of the three Köppen-Geiger climate zones, Cfb, Csa, and Dfb) adaptation to high temperatures of the population sampled in Israel, consistent with the results obtained for another Israeli population investigated by Zhan and McDonald (2011);and (iii) (Suffert et al., 2015). This study thus reveals a two-tier thermal adaptation, with seasonal dynamics nested within and potentially occurring in each geographic local adaptation over annual epidemics. This key finding shows that adaptive patterns are "eco-evolutionary snapshots" that should be interpreted with caution, to such an extent that certain evolutionary dynamics of microbial populations can be of one type over a very short timescale and another type over longer timescales. Indeed, adaptive dynamics may differ with the timescale investigated (annual or pluriannual), particularly for annual crop pathogens with both sexual and asexual reproduction cycles, such as Z. tritici . Our findings could be summarized by the counterintuitive statement "local seasonal adaptation is stronger but more fleeting than geographic adaptation" although we would expect that regions with lower seasonal contrasts in temperature (e.g., with mild winters) will exert weaker selective pressure. The use of sequential temporal sampling would make it possible to capture shifts in thermal adaptation over and between wheat-growing seasons and to detect potential trade-offs between aggressiveness and survival over winter (e.g., Montarry et al., 2007).

| From adaptation patterns to ecoevolutionary processes
Consistent with previous studies, our findings highlight the existence of high levels of genetic diversity and an absence of its structuring across Z. tritici populations collected from local wheat fields (Zhan et al., 2001) up to the regional and continental scales (Linde et al., 2002;Schnieder et al., 2001) or over the course of an epidemic cycle (Chen et al., 1994;Morais et al., 2019). The high level Categorization tackling functional redundancy at the population level (i.e., whether the thermotypes composing a population are more or less well differentiated within the whole functional space). The three populations presented here demonstrate the relevance of considering functional redundancy vs. vacant functional space when assessing emergent properties of populations such as a generalist nature at population level (e.g., a generalist population can be composed of specialist individuals with narrow individual TPB 80 distributed over the functional space, resulting in broad TPB 80 population coverage). (d) The translation of populations into functional groups makes it possible to investigate group-level selection, testing for general assumptions of adaptation to given environments (e.g., competitive advantage of lowtemperature specialists in cold environments, generalists in variable environments, and high-temperature specialists in warm environments) providing insight into the potential of populations to adapt to changes in their environment (a subtle balance between diversity levels for intra-and interindividual variation in thermal responses) high diversity of thermal responses in Z. tritici highlighted here, their heritability (Lendenmann et al., 2016), and the high level of local heterogeneity within wheat canopies (Chelle, 2005) suggests that local thermal conditions probably exert strong selection pressure on thermal sensitivity (for which TPCs are probably the best proxy as they may themselves be direct targets of selection; Scheiner, 1993;Via, 1993), even in the presence of high gene flow due to long-distance ascospore migration (

| Functional group composition: an operational approach for investigating population dynamics
Our study illustrates how the functional classification of TPCs into thermotypes with multivariate statistical procedures can provide a complementary means of deciphering diversity patterns in the biological responses quantified in reaction norms. In particular, it constitutes an operational tool for assessing functional similarity at the individual level (i.e., whether the apparent variation observed in thermal parameters is functionally significant; Figure 8a,b) and at the population level (i.e., whether the thermotypes constituting a population are more or less well differentiated within the whole functional space; Figure 8c). However, caution will be required in the extension of this approach to comparisons over multiple data sets, through the development of comparable classification systems, taking into account the variation of the classification with the populations sampled by explicitly stating which ranges of trait values are hidden behind a given group description (e.g., affixing levels of "adaptation" in the sense of higher performance at given temperature ranges: very low, low, high, very high). This description of populations in terms of functional groups makes it possible to move from a description of phenotypic patterns and shifts in population composition to an inference process. This process may, for example, be based on comparisons of the competitive advantage of thermotypes under given thermal scenarios: for example, "do more variable environments favor thermal generalists?" or "is there a shift in the optimal range of thermal responses with mean temperature conditions?" (Figure 8d). This classification into thermotypes enabled here to go beyond a purely descriptive framework, and future investigations will need to be undertaken to tackle the physiological bases of these differentiations in thermal responses. The thought-provoking results of Francisco et al. (2019) could be used to test whether several strains belonging to those thermotypes also correspond to specific or main morphotypes that would increase their tolerance under some environmental conditions (e.g., if "warm-adapted" individuals exhibit higher proportions of stress-tolerant growth forms such as chlamydospores under warmer temperatures). All in all, this functional approach lays the foundations for future studies of the potential of pathogen populations to adapt to changes in their environment, from seasonal changes in the short term, to global warming in the long term. In particular, it will prove useful in gaining a fuller understanding of how new aggressive fungal strains may emerge and expand into previously unfavorable environments (Milus et al., 2009;Mboup et al., 2012;Stefansson et al., 2013). This is a crucial area of investigation that is all too often overlooked in models for predicting plant disease epidemics in conditions of climate change (West et al., 2012).

| CON CLUDING REMARK S
The detailed characterization of a microbial phenotype as a profile rather than a mean allowed for analyses that accounted for the range of sensitivities of individual strains rather than solely their mean sensitivity. This gave insight into the high level of functional divergence in the plasticity and variation of individual thermal responses over geographic and seasonal scales, highlighting the occurrence of two-tier dynamics in thermal adaptation. These findings raise intriguing questions regarding the mode of selection operating on these functional groups of individuals with similar competitive advantages in given thermal conditions. Deciphering the mechanisms underlying this maintenance of diversity in population phenotypic composition will prove useful for expanding our understanding of eco-evolutionary responses and the potential of populations, species, and communities to adapt to environmental change.

ACK N OWLED G M ENTS
We would like to thank Aigul Akhmetova ( leaves with STB symptoms, which was required for the establishment of the Euro-Mediterranean Z. tritici collection upon which this study is based. We would also like to thank Sylvain Pincebourde (IRBI-CNRS, France) for many fruitful discussions.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The thermal phenotyping data set associated with this manuscript is