Mapping the adaptive landscape of a major agricultural pathogen reveals evolutionary constraints across heterogeneous environments

Pathogen success in heterogeneous biotic and abiotic environments depends on genetic variation in life-history traits that govern host infection, reproduction, dispersal, and survival. Adaptive evolution can be constrained by trade-offs arising from differential resource allocation and antagonistic gene actions (including pleiotropy) that prevent simultaneous increases in multiple favorable traits. Genetic correlations among traits may be generated by pleiotropy or by linkage disequilibrium, with pleiotropy expected to impose more stable evolutionary constraints. Most fitness-related traits are polygenic, complicating the identification of trade-offs and pleiotropic effects. Zymoseptoria tritici, the causal agent of septoria tritici blotch on wheat, has a complex life cycle with adaptations for growth across temperatures, asexual reproduction via pycnidia, and survival structures for off-season persistence. Prior phenotypic work suggested correlations and potential trade-offs between lesion development (virulence) and pycnidia production (reproduction). Global populations show variation in thermal adaptation and fungicide resistance, and the pathogen exhibits variable melanin production under stress and chlamydospore formation at high temperature. However, the genetic architecture of key life-history traits and how it gives rise to trade-offs and pleiotropy remain largely unknown. The study aims to map the genetic basis of 50 life-history traits across diverse host and non-host environments in a global collection of Z. tritici, quantify heritability, and resolve genetic correlations to reveal constraints and facilitations shaping pathogen evolution.

The study builds on theory and evidence that trade-offs and pleiotropy constrain adaptation. Prior work demonstrated environment-dependent trade-offs: thermal trade-offs in Escherichia coli, fitness costs of drug resistance in Candida albicans, and a link between melanin synthesis and virulence in Aspergillus fumigatus. In Z. tritici, phenotypic correlations have been noted between lesion development and pycnidia production, melanin production and growth/fungicide susceptibility, as well as the formation of chlamydospores under high temperature. Despite these insights, mutations and genomic loci underlying synergistic or antagonistic interactions among traits (virulence, reproduction, stress tolerance) have largely remained unidentified, and most fitness-related traits are expected to be polygenic, necessitating genome-wide approaches to dissect trade-offs and pleiotropy across many traits.

Sampling and isolates: A global panel of 145 genetically distinct Zymoseptoria tritici isolates was assembled from five field populations across four continents (Australia n=27; Israel n=30; Switzerland two nearby sites n=2 and n=30; USA Oregon populations sampled from two cultivars: Madsen [Oregon.R] n=26 and Stephens [Oregon.S] n=30). Collections span 1990–2001. Isolates were stored at -80 °C. Whole-genome sequencing and SNP calling: Illumina HiSeq2000 paired-end (100 bp, ~500 bp insert) data were available for 130 isolates and generated de novo for 15 additional isolates (BioProject PRJNA327615). Reads were quality-trimmed with Trimmomatic v0.36, aligned to the IPO323 reference genome using Bowtie2 v2.3.3, and duplicates removed with Picard MarkDuplicates. GATK v4.0.1.2 HaplotypeCaller (haploid, GVCF mode) followed by joint genotyping (GenotypeGVCFs, maxAltAlleles 2) was used. Variants were filtered with stringent thresholds: QUAL<250, QD<20.0, MQ<30.0, BaseQRankSum not in [-2,2], MQRankSum not in [-2,2], ReadPosRankSum not in [-2,2], FS>0.1. Final dataset: 716,619 high-quality biallelic SNPs (call rate ≥80%, MAF>5%) across 21 chromosomes. In planta phenotyping: Virulence (necrotic lesion area) and asexual reproduction (pycnidia density within lesions) were quantified on 12 wheat genotypes (five landraces: Chinese Spring, 1011, 1204, 4391, 5254; six cultivars: Drifter, Gene, Greina, Runal, Titlis, Toronit; one backcross line ArinaLr34). Greenhouse conditions: 22/18 °C day/night, 70% RH, 16-h photoperiod. For each isolate, uniform spray inoculation to runoff was performed; three independent biological replicates in different chambers across two experimental phases. Leaves were harvested 19–26 dpi and analyzed by automated image analysis (AIA) to derive lesion area and pycnidia density. Host specialization index for reproduction per isolate was computed as the adjusted coefficient of variation of LSmeans across hosts. In vitro phenotyping: Growth and stress responses were measured on PDA at 15 and 22 °C from mean colony area across time (8, 11, 14 dpi) using AIA; growth rates (mm/day) were estimated by linear models on mean colony radius over time. Temperature sensitivity was calculated as growth rate ratios (15 vs 22 °C). Fungicide sensitivity was assessed via growth at 22 °C on PDA with propiconazole (0.05 ppm) and via microtiter dose-response assays (propiconazole series from 0.00006 to 1.5 mg/L) to estimate EC50 using the ‘drc’ package. Melanization was quantified from grayscale values (converted to 0–255 black scale) at 8, 11, 14, 18 dpi under 15 °C, 22 °C, and 22 °C with fungicide (0.05 ppm). Morphological stress response (chlamydospores/hyphae formation) followed Francisco et al. Data were based on multiple technical replicates (e.g., ~45 colonies per isolate for growth). Statistical analyses: For in planta traits, log-transformed LSmeans per host×isolate from Dutta et al. were used. For in vitro traits, LSmeans per isolate were extracted via one-way ANOVA. A clustered heatmap of z-scored trait values was generated using ComplexHeatmap. GWAS was conducted for 50 traits using 716,619 SNPs across 145 isolates. Population structure was explored via PCA; a genomic relatedness matrix (GRM) was computed using normalized identity-by-descent in TASSEL. Based on BIC, principal components were not included as fixed covariates; GWAS used a mixed linear model (MLM+K, GRM as random effect) in GAPIT (model: y = Xb + Zu + e). Multiple-testing control used Bonferroni (α=0.05) and FDR (5% and 10%; Benjamini-Hochberg). SNPs surpassing FDR 10% were annotated by nearest-gene mapping (BEDtools ‘closest’) against the IPO323 annotation; protein features were assigned by InterProScan; secretion signals predicted with SignalP, Phobius, TMHMM; GO enrichment tested with GOstats (FDR 0.01, minimum term size 5). Heritability and genetic correlations: SNP-based narrow-sense heritability (h2_snp) was estimated per trait using GCTA GREML with the GRM as random effect (haploid adaptation noted). Genetic correlations (rg) between all trait pairs were estimated via bivariate GREML, deriving genetic covariance and variances from the GRM. Genetic correlation networks were visualized with qgraph. Phenotypic correlations used Pearson correlations on z-scores with BH correction; binary morphological response correlations used point biserial correlation.

- Extensive quantitative variation in 50 life-history and stress-response traits was observed across 145 Z. tritici isolates from four continents, with substantial within-population variability and differing mean trait values among populations. - Population-level patterns: The Israeli population showed the highest host specialization for reproduction and the highest proportion of isolates forming chlamydospores under stress; the Swiss population had the highest azole fungicide resistance, consistent with earlier and more intense azole use in Europe. - Most traits exhibited polygenic architectures with high SNP-based heritability, indicating that short-term adaptation is likely to proceed via allele frequency shifts at many loci rather than single large-effect mutations. - Heritability estimates: Virulence h2 ranged from 0 to 0.59 (SE ≈ 0.14), whereas reproduction h2 ranged from 0.43 (SE ≈ 0.15) to 0.91 (SE ≈ 0.03). Mean h2 for reproduction was higher (mean ≈ 0.72, SE ≈ 0.10) than for virulence (mean ≈ 0.47, SE ≈ 0.15), suggesting stronger and faster potential response to selection for reproduction. - Genetic correlation analysis revealed negative genetic correlations between traits involved in host colonization/damage (e.g., virulence) and survival or performance under stressful abiotic conditions, indicating antagonistic pleiotropy and evolutionary constraints on maximizing both simultaneously. - Conversely, adaptation to certain abiotic stressors (e.g., temperature extremes, fungicide exposure) appeared to be aided by synergistic pleiotropy among stress-response traits, potentially facilitating coordinated improvement under stress. - No strong geographic structuring of overall phenotypic profiles was detected, despite population differences in specific trait means. - The integrative GWAS and correlation framework mapped adaptive loci across the genome and connected them to trait interdependencies, supporting predictions of evolutionary trajectories under environmental perturbations.

The study addresses how genetic architecture and pleiotropy shape the adaptive landscape of a major crop pathogen across diverse environments. High heritability estimates, particularly for reproduction, imply that selection on reproductive capacity can produce rapid evolutionary change. The prevalence of polygenic control suggests adaptation will largely involve coordinated allele frequency shifts across many loci. Negative genetic correlations between host damage traits and stress-survival traits provide evidence for antagonistic pleiotropy, imposing constraints that may prevent simultaneous optimization of virulence and stress tolerance. Such constraints could limit epidemic potential under certain environmental regimes or management interventions. In contrast, positive genetic correlations among abiotic stress-response traits indicate synergistic pleiotropy that may facilitate adaptation to temperature and fungicide pressures. Together, these findings help explain observed population differences (e.g., higher fungicide resistance in Europe) and offer a genome-informed basis for predicting pathogen evolutionary responses to climate variability and control measures.

By combining large-scale GWAS with comprehensive phenotyping across host and non-host environments, the study maps the polygenic architectures and genetic correlations underlying 50 life-history traits in Z. tritici. It reveals strong heritability for many traits, especially reproduction, and identifies antagonistic pleiotropy between virulence-related traits and stress survival, alongside synergistic pleiotropy among abiotic stress responses. These insights clarify evolutionary constraints and facilitations that shape pathogen adaptation and can inform risk assessments and management strategies under environmental change and fungicide use. Future work could functionally validate candidate loci, assess temporal dynamics of allele frequency shifts in field populations, expand host panels and environmental conditions, and integrate gene-by-environment interactions to refine predictive models of pathogen evolution.

- Sample size, while substantial for a pathogen GWAS (n=145), was acknowledged as relatively small, potentially limiting power to detect small-effect loci and complex interactions. - Phenotyping was conducted under controlled greenhouse and laboratory conditions, which may not capture the full spectrum of field environmental variability and host–pathogen interactions. - Genetic correlations estimated via GREML rely on the GRM as an approximation of causal variant relationships; accuracy depends on SNP density and linkage to causal sites. - Some analyses excluded binary traits (e.g., morphological stress response) from certain multivariate visualizations, potentially reducing integrative comparisons. - Reliance on a single reference genome (IPO323) and current annotations may miss accessory genomic variation or structural variants contributing to traits.