A thousand-genome panel retraces the global spread and adaptation of a major fungal crop pathogen

Global genetic structure of the pathogen tracks the historical spread of wheat

We assessed the evolutionary trajectory of the pathogen in conjunction with the history of global wheat cultivation (Fig. 1a). For this, we assembled a worldwide collection of Z. tritici isolates from naturally infected fields (Fig. 1b). We selected isolates covering most wheat production areas, both in the center of origin of the crop (i.e., the Fertile Crescent in the Middle East), and in areas where wheat was introduced during the last millennia (i.e., Europe and North Africa), or last centuries (i.e., the Americas and Oceania; Fig. 1c). We called variants in a set of 1109 high-quality short-read resequencing datasets (Supplementary Data 1, 2) covering 42 countries and a broad range of climates. Using a joint genotyping approach, we produced raw variant calls mapped to the telomere-to-telomere assembled reference genome IPO323. To assess genotyping accuracy, we used eight isolates with replicate sequencing data to analyze discrepancies. We adjusted quality thresholds targeting specifically the type of genotyping errors observed in our data set (Fig. S1). The improved filtering yielded 8,406,818 high-confidence short variants (short indels and SNPs). The final variant set included 5,578,488 biallelic SNPs corresponding to 14.1% of the genome.

We tested whether global diversity patterns of pathogen populations are likely a consequence of the history of wheat cultivation. We first performed unsupervised clustering of genotypes and identified eleven well-supported clusters (Fig. 2a, Figs. S2,3). Over 90% of the genotypes were clearly assigned to a single cluster (Fig. 2a, Supplementary Data 3). Two clusters were identified among genotypes originating from the pathogen center of origin, distinguishing collections from Iran and Middle Eastern regions. Genotypes from Africa and Europe split into two distinct genetic clusters without any apparent secondary structure within clusters. This lack of any fine-scale structure is remarkable given the extensive geographic sampling of European genotypes and suggests extensive gene flow within the continent. Genotypes from Oceania grouped into three distinct clusters marked by collections from Tasmania, the Australian mainland, and New Zealand. Genotypes from North America formed two clusters along a North-South separation. Finally, South American genotypes formed two clusters split along the Andes (Chile versus Argentina and Uruguay). Some uncertainty exists in the assessment of regional population structure by low coverage of major wheat-producing countries such as Russia and Ukraine. Septoria tritici blotch is only sporadically reported in China. In complementary analyses, we found that a phylogenetic network accounting for the high frequency of recombination consistently reflected the global population structure (Fig. S4). A principal component analysis of all genotypes confirmed the nested genetic structure with differentiation at the continent level, subdivisions within some continents and the existence of admixed genotypes (Fig. 2b, Fig. S5).

We analyzed the history of population splits and admixture using allele frequency information (Fig. 2c). The analyses largely supported a genetic structure shaped by the introduction of wheat across continents. The historical relationships between clusters show an early divergence of the Middle Eastern and North African clusters matching the early introduction of agriculture in these regions. Populations in Europe and the Americas share a similar time point of divergence consistent with extensive contributions of European genotypes to the Western hemisphere. Oceanian groups have diverged as a single branch from genotypes most closely related to extant European populations. Matching the introduction of wheat to Oceania from the European continent, the Australian and New Zealand pathogen populations share a common origin rooted in European genetic diversity. Populations from Australia show also a striking loss of diversity and higher linkage disequilibrium compared to European diversity consistent with a significant founder effect (Fig. 2d, e). Similarly, populations in South and North America have reduced genetic diversity compared to extant European populations as suggested previously based on Sanger sequencing¹⁶. The highest diversity was found in populations from Africa and the Middle East closest to the center of origin. Overall, the global genetic structure of the pathogen reveals multiple founder events associated with the introduction of wheat to new continents.

Ongoing gene flow among regions should lead to admixed genotypes. We found that nearly 10% of all analyzed genotypes showed contributions from at least two clusters. The most significant recent gene flow was detected between Middle Eastern/North African clusters and European clusters in North Africa (i.e., Algeria and Tunisia) as well as Southern and Eastern Europe (i.e., France, Italy, Hungary, Ukraine, Portugal, and Spain; Supplementary Data 3). We found a particularly high incidence of recent immigration in a durum wheat population in the south of France. The population consisted only of hybrids or atypical genotypes suggesting either recent migration from North Africa or host specialization on durum wheat varieties. Additionally, we found hybrid genotypes with European ancestry in both North America and in Oceania. The relatively balanced ancestry proportions in these hybrids suggest very recent gene flow dating back to only a few generations. We further investigated past gene flow between clusters by allowing Treemix to infer migration events, thus creating a population network (Fig. S6a–d). Three distinct recent migration events were best explaining the data with specific migration routes from the Middle East/African clusters to North America, from an Australian cluster to South America and between two Oceanian clusters (Fig. S6d). However, the migration events did not affect the overall shape of the inferred population tree (Fig. 2c, Fig. S6b–d). To better understand effects of long-distance gene flow, we investigated the relationship between relatedness among genotypes (i.e., identity-by-state) and geographic distance. At the continent level, we observed a negative relationship between identity-by-state and geographic distance (Fig. S7). The wide distribution of identity-by-state values shows that although closely related isolates tend to be found at closer geographic distance, distantly related isolates can be found at both far and close geographic distances. Long-distance migration events are most likely caused by international trade similar as for other crop pathogens^17,18,19. In combination, our findings show an important role of long-distance dispersal impacting the genetic make-up of populations from individual fields to continental scale genetic diversity.

Relaxation of genomic defenses against transposable elements concurrent with global spread

Transposable elements (TEs) are drivers of genome evolution. In Z. tritici, TE activity created beneficial mutations for fungicide resistance and virulence on the wheat host^20,21. Rapid recent adaptation of the pathogen has benefitted from the activity of TEs with consequences for genome size²². Unchecked transposition of TEs can be deleterious and an array of defenses mechanisms has evolved to counteract their activity both at the genomic and epigenetic level including targeted mutations and silencing²³. To analyze the effectiveness of genomic defenses against active TEs, we screened all genomes for evidence of TE insertions. We mapped short-read sequencing data on the reference genome and a species-specific TE sequence library. We classified evidence for TEs in each of the analyzed isolates as reference TEs (i.e. also present in the reference genome) and non-reference TE (i.e. absent). Detected TEs among isolates were binned into loci (width 100 bp) to account for uncertainties about the precise mapping of the insertion point. We found that the frequency spectrum of TE insertions is heavily skewed towards low frequencies with 77% of TE insertions being found in single isolates (~0.1% frequency) and 96% of insertions were found in ten or fewer isolates (<1%; Fig. S8a) consistent with strong purifying selection. The Z. tritici genome contains both core and accessory chromosomes (i.e., chromosomes not shared among all isolates of the species; Fig. 3a)²⁴. Accessory chromosomes have higher TE densities (15–33% vs. 5.5–19%²⁴; Fig. 3b) reflecting lower selection pressure on accessory chromosomes²⁵. Beyond this, accessory chromosomal regions are broadly differentiated from core regions based on sequence, transcription and epigenetic feature sets (Fig. 3c). The primary differentiator (i.e., the first principal component) separated euchromatic and heterochromatic (H3K9me2) chromosomal segments. These differences in epigenetic marks were matching TE density variation, consequences of genomic defenses (i.e., repeat-induced point mutations), and GC content. The number of TE insertions was not correlated with GC content but was positively correlated with the presence of facultative heterochromatin (H3K27me3) and negatively correlated with euchromatin marks. H3K27me3 is also a hallmark of accessory chromosome segments consistent with previous findings^26,27,28. Short insertion/deletion (indel) polymorphism was positively correlated with the heterochromatic accessory regions, consistent with purifying selection acting against indels mostly in gene-dense regions.

The pathogen genome shows variation in TE content among populations from different continents, suggesting that the TE content was influenced by the colonization history of the species^22,29. Hence, we tested whether recently established populations were differentiated from center-of-origin populations. We found that TE insertion polymorphisms (TIPs) are most often specific to a genetic cluster with only 31% of insertions shared by two or more populations (Fig. S8b). We used TIPs as a genetic marker and found that population differentiation was matching differentiation assessed using short variants (Fig. S8c). Hence, the population history of the pathogen was an important factor shaping TE content. We also found that the per-individual TE content has increased in most areas outside of the center of origin. TE detection using short read sequencing data is constrained by sequencing depth, hence we accounted for depth in our comparative analyses of TE content (Fig. 3d). The TE content of the Middle Eastern and African genomes was lower than in any other regions ranging from 164–394 (median = 302) and 184–429 TIPs (median = 300), respectively, including only samples with coverage below 35. In contrast, the Oceanian, North American and South American genomes contained among the highest TE numbers (258–513, median = 373; 283–513, median = 462; 237–476, median = 375 TIPs, respectively) and European genomes showing intermediate TE counts (209–669, median = 334 TIPs, including only samples with a coverage below 35).

The increase in TE content was broadly spread among superfamilies and included both DNA transposons and retrotransposons (Fig. S9, 10). Some superfamilies showed particularly strong signatures of expansion outside of the Middle East and Africa including Copia retrotransposons (RLC) and other LTR retrotransposons (RLX) (Fig. 3e, f). The RLC expansion is mainly driven by the Deimos element with an increase from 11.3 copies on average in genomes from the Middle East and Africa to an average of 50.8 copies per genome outside the center of origin (Fig. 3g and S9, 10). The RLX superfamily expansion is driven by the Styx element, generally more present outside the center of origin (Fig. 3g and S9,10). The Styx activity was shown to affect asexual reproduction on the host³⁰. We also observed population-specific expansions such as the LINE (RIX) Juno element mostly present in North American and Australian populations, the Helitron (DHH) Fatima expanded primarily in the Australian population, and the DTX miniature inverted-repeat TE (MITE) Unicorn most widespread in North American and New Zealand populations (Fig. S9,10). The broad increase in TE content outside of the center of origin suggests a relaxation of genomic defenses over the evolutionary history of the pathogen on wheat.

Many ascomycetes share a genome defense mechanism against TEs that can rapidly introduce targeted mutations into newly duplicated sequences, called repeat-induced point mutation (RIP)^31,32. RIP machinery is active in genomes of Z. tritici, with high levels of RIP-like mutations identified in genomes from the center of origin and wild-grass-infecting sister species^29,33. We analyzed the global panel for evidence of RIP-like mutations by reporting the RIP composite index. The median index is above 1 in TE sequences across all genetic clusters. However, we found that RIP strength varies considerably among genetic clusters with the strongest signatures found in genomes from Middle Eastern and African isolates (Fig. 4a). The Middle Eastern and African clusters tend to include genomes with both low TE content and strong RIP signatures (Fig. 4b, S8a, b). In more recently colonized regions, genomes showed a negative correlation between the strength of RIP signatures and the amount of TEs per genome (p value <2.2e−16; slope of 6.7 × 10⁻⁴ in all populations outside of the Middle-East and slope of 4.6 × 10⁻⁴ in the European population only; Fig. 4b & S11a, b). The association between TE content and strength of RIP signature is consistent with genomic defenses being less capable to prevent new TE insertions following migration out of the center of origin.

How genomic defenses against TEs are modulated in fungi remains largely unknown. The RIP machinery is activated during sexual recombination in Neurospora³⁴, thus suggesting that reduced sexual reproduction could lower the efficacy of RIP. In Z. tritici, sexual reproduction can occur during the wheat growing season but most reproduction is thought to occur at the end of the season³⁵. The rate of sexual reproduction is high in all populations as we found that the ratio of mating types was consistently close to 50:50 (Fig. S11c). Alternatively, the RIP machinery could have lost function in some populations. In this case, TEs should show bimodal signatures of RIP with old TEs carrying signatures of historical RIP activity and recent insertions with no or weak RIP signatures. Indeed, a high percentage of TEs shows only weak evidence for RIP (composite index <0.5) in genetic clusters outside of the Middle East and Africa (mean = 20.1%, range from 11.3–34.6%; Fig. S11d). By contrast, in the center of origin all genomes had less than 10% of TEs showing weak evidence for RIP (mean = 7.48%, min = 5.18%). We confirmed the bimodal RIP signatures by analyzing individual TEs in 20 chromosome-level assembled genomes (Fig. 4c). All genomes shared a major RIP composite index peak of ~2. A secondary peak indicative of TEs without RIP was found in genomes from the Americas, Oceania, and Europe. Despite the higher TE activity, we found no evidence that the loss in RIP functionality also led to high gene duplication events. The strength of RIP signatures is associated with the size of the TEs as small elements show little to no RIP signatures (Fig. 4d) consistent with the RIP machinery being inactive on small repeats²³. In Neurospora crassa, two pathways mediate RIP with one dependent on RID and one on Dim2³⁴. In Z. tritici, the Dim2 methyltransferase was functionally linked to the occurrence of RIP-like mutations in repeats including during mitosis^33,36. Furthermore, the presence of a functional dim2 copy is strongly correlated with lower GC content of TEs, hence deamination of cytosines^33,36. The dim2 gene was duplicated multiple times in some genomes causing loss-of-function mutations triggered by the RIP mechanism itself^33,37. We found that the ancestral copy of dim2 shared higher identity with the functional copy of dim2 in the genomes of the Middle Eastern isolates than other populations (Fig. S11e). The European populations had both the largest range in the number of detected paralogs and the highest copy numbers overall (Fig. S11f). This is consistent with a deleterious runaway gene-duplication process affecting a molecular component of the RIP machinery, explaining the loss of RIP efficacy within the species.

Adaptation to fungicides and changing climates along continental gradients

Globally distributed pathogens experience significant environmental heterogeneity that potentially constrains or facilitates future range expansions and adaptation. The use of pesticides across the globe to combat agricultural pathogens has triggered the parallel emergence of resistance with significant economic consequences⁶. To retrace the global emergence of fungicide resistance in Z. tritici, we analyzed mutations in resistance genes using isolates collected over three decades (1986 to 2016). This time span covers the introduction of several major fungicides to agricultural fields. Resistance to the ubiquitously used azole fungicides is often mediated by mutations in the CYP51 gene³⁸. Recent North American populations gained the Y137F mutation but not the I381V or V136A mutations rising in frequency in Europe over the same period (Fig. 5a). European populations harbored the most diverse set of azole resistance mutations, consistent with the early and intense applications of this fungicide class. The I381V and V136A mutations occurred at high frequency since the early 2000s, whereas the S524T mutation was only observed at a low frequency with a delayed onset. Resistance arose later in Oceania and North America, consistent with the later application of azoles in those locations. No resistance mutations were detected in the Middle Eastern or African populations, matching the absence of azole treatments in those regions. We found similar geographic patterns for the E198AK mutation in the beta-tubulin gene associated with benzimidazole resistance as well as for the G143A mutation in the mitochondrial gene cytb, known to cause resistance to Quinone outside inhibitors fungicides (Fig. 5b). As expected from their more recent introduction, mutations related to resistance to succinate dehydrogenase inhibitors (SDHI) were only observed in the most recently sampled populations in Europe (Fig. S12). Overall, the global analyses of fungicide resistance signatures show how European populations consistently developed the first known mutations to newly introduced fungicide.

Changes in climatic conditions create complex challenges for plant breeders to create resilient crops³⁹. Concurrently, pathogen populations are exposed to changes in temperature and humidity patterns. The historic spread of Z. tritici has likely created significant selection pressure to adapt to climates associated with the global range of wheat cultivation. Here, we analyze the genetic architecture of climate adaptation by mapping standing variation along climatic clines. The pathogen is endemic to regions with distinct climates, from the dry and warm conditions in the Middle East to the temperate oceanic climate of Western Europe and the humid continental climates of some North American locations. We performed genotype-environment association (GEA) mapping based on the climatic conditions of the sampling locations. We analyzed a total of 19 bioclimatic variables covering annual trends, seasonality, and extreme environmental factors, such as the maximum temperature of the warmest month or the precipitation of the driest quarter year (Fig. 5c, d, Supplementary Data 4).

We identified 1956 variants significantly associated with at least one climatic variable and 640 variants with a minor allele frequency (MAF) higher than 5%. The number of associated variants per climatic variable ranged between 36 and 541 (for BIO9 and BIO6 respectively), including 1–190 significant SNPs with a MAF > = 5% per climatic variable (Fig. 5e). We investigated whether variant classes were enriched in the set of significantly associated SNPs. Using permutations, we found that both non-synonymous and intergenic variants were more frequently associated with climatic variables than expected randomly while synonymous variants were significantly depleted (Fig. S13b). We found 187 genes that were in proximity or directly affected by variants associated with bioclimatic variables and with a MAF > = 5%, including 65 containing non-synonymous variants (Supplementary Data 5). For each GEA, we retrieved significantly associated loci by clustering significant variants within a distance of 10 kb. The significant variants clustered into 5–27 distinct loci per GEA consistent with a polygenic basis for most climate adaptation (Fig. S14–19; Supplementary Data 4). A polygenic architecture of thermotolerance was previously found in other fungi such as Saccharomyces cerevisiae^40,41. A large number of associated loci were shared between GEA of different climatic variables. Highly correlated climatic variables including BIO5 and BIO10, the maximum temperature of the warmest month and the mean temperature of the warmest quarter, respectively, shared also higher numbers of significantly associated SNPs (Fig. 5d, Fig, S13a). However, we also identified hotspots of climatic adaptation loci for largely independent climatic variables (Fig. 5f). We identified a large segment of chromosome 7 and a telomeric region of chromosome 13 to be hotspots for climate associations. The chromosome 7 locus overlaps with the Cyp51 gene involved in azole fungicide resistance. Hence, the association could be due to correlations in the application of fungicides and climatic factors. Temperate regions such as Europe show higher azole resistance (see above) than the Middle East, thus leading to an indirect association between fungicide resistance genes and climatic variables. However, at the same location on chromosome 7 is a quantitative trait locus (QTL) for growth at suboptimal temperature⁴², so that locus could well underpin climatic adaptation independent of Cyp51. Further analyses of temperature sensitivity QTLs showed that three out of the four previously described loci are overlapping with loci associated with climatic variables (Supplementary Data 6). For example, the temperature-sensitivity QTL on chromosome 1 overlaps with loci associated with multiple temperature-related climatic variables such as the mean temperature of the warmest quarter (Fig. 6a–c). A variant associated with the mean temperature of the warmest quarter and overlapping with a QTL previously discovered in a central European cross (chr 1 at 2,090,068 bp) shows a global distribution of both alleles. By contrast, a second variant (chr 10 at 452,864 bp) associated with the mean temperature of the warmest quarter but not overlapping the previously discovered QTL (Fig. 6b, d), shows no allelic variation within central European populations. Alleles associated with warmer climates were also found in the US Midwest characterized by cold winters and high temperatures in summer. Our association analyses capture most likely only a subset of all loci contributing to climatic adaptation across the global distribution range. Hence, some degree of climate-genotype mismatches is expected across geographic regions if climatic adaptation is polygenic. Our analyses highlight the power of covering global genetic diversity to gain insights into the genetic architecture of recent adaptation in species.

To identify whether adaptive mutations tend to arise locally or occur across large geographic areas, we clustered significant loci based on the presence or absence of adaptive variants per country. We identified eight clusters characterized by shared adaptive variants with a similar geographic distribution (Fig. 6e, Fig. S19). Some clusters of adaptive variants are highly geographically localized (i.e., clusters 3, 5, and 6) while other clusters are widespread (i.e., cluster 4). Most adaptive alleles for extreme cold conditions were found in the populations subjected to the harsh winters of continental North America (see cluster 3; Fig. 6e–g). Clusters of adaptive alleles were geographically widespread such as the distributions along the Mediterranean coast (see cluster 1). Taken together, the global genome panel revealed substantial standing variation for environmental adaptation with complex geographic patterns of local adaptation evolution.

Discussion

Analysis of a thousand-genome panel recapitulated the spread of a major fungal crop pathogen revealing tight links to the history of global wheat cultivation. The early divergence between Middle Eastern and African genotypes from those collected in the rest of the world is consistent with a single expansion event from the center of origin dating back millennia. The extant genetic variation was strongly shaped by successive colonization bottlenecks during the introduction of the pathogen to the Americas and Oceania. The distinct loss of genetic diversity and increased linkage disequilibrium likely caused the loss of adaptive genetic variation and reduced evolutionary potential in the most recently colonized regions. TE activity has underpinned the rise of major adaptive mutations in the species^20,30,43. Remarkably, the TE activity is underpinned by a marked relaxation or even loss of genomic defenses following population bottlenecks during global colonization. The higher activity of TEs is likely a direct consequence of reduced control and may have long-term consequences for the pathogen. The relaxation of genomic defenses in populations from the American, Oceanian, and European continents could have been selected as an evolvability trait, thus increasing variance in fitness in populations. The relaxation of genomic defenses likely underpins incipient genome expansion within the species while increasing the risk of mutational meltdown. The resilience of crops and agricultural ecosystems is threatened by the changing climate. The ability of pathogens to adapt and expand their range under altered humidity and temperature regimes as well as changes in seasonal patterns is a major concern. The identified genomic regions associated with adaptation to environmental conditions highlight how a global pathogen carries extensive variation to cope with climate change. Integrating population genetic information of pathogen adaptation to climatic gradients is a powerful asset for risk models of future pathogen spread. As climate change impacts crop production and most likely resistance to pathogens, future research explicitly addressing joint host and pathogen responses to warmer and drier climates is necessary to mitigate future crop losses.

Source: Ecology - nature.com