Genome-wide Evidence of Host Specialization in Wild and Farmland Populations of the Fungal Leaf Spot Pathogen, Cercospora beticola

Abstract

One of the most recent crop species to be domesticated is sugar beet (Beta vulgaris L. ssp. vulgaris Doell.), which was bred for high sucrose content within the last few centuries in Europe. Crop domestication can also lead to the evolution of novel pathogens, which may spread across large geographical distances with their crop host. In this study, we addressed the recent evolution of the fungal pathogen causing the disease Cercospora leaf spot, Cercospora beticola. This pathogen has become increasingly important in sugar beet and table beet production worldwide. We used genome sequences of 326 C. beticola isolates collected from 4 continents from 4 closely related Beta subspecies (3 domesticated and 1 wild). We applied population genomic analyses to identify signatures of population differentiation and host specialization in C. beticola populations derived from the cultivated and wild hosts. We found evidence that C. beticola populations in agro-ecosystems likely originate from sea beet-infecting isolates. Intriguingly, host jumps from wild to cultivated beet occurred in at least 2 independent events as evidenced by our population data of C. beticola from wild beet collected in the Mediterranean and the UK. We explore the occurrence of genetic variants associated with fungicide resistance and virulence and show that standing genetic variation in C. beticola populations from both wild and domesticated plants may serve as a reservoir of functionally important alleles. Overall, our results highlight the ability of C. beticola to invade the agro-ecosystem and establish new populations, demonstrating the rapid adaptation potential of the species.

Introduction

Population genetics has proven to be a valuable tool in our efforts to understand the evolutionary trajectories of economically significant plant pathogens at a global scale (Pereira et al. 2020; Thierry et al. 2022; Feurtey et al. 2023; Taliadoros and Stukenbrock 2023). These studies provide valuable insights into the adaptation of pathogen populations to local conditions and hosts, which, in turn, offer crucial information for disease management, measures for regulating the movement of plant materials, and the use of fungicides to minimize disease outbreaks. For instance, insights into the origin of pathogen species, patterns of dispersal, processes of host specialization, and mechanisms of evolution and dispersal have been achieved by sampling of pathogen populations from different geographical regions and different hosts and by detailed analyses of polymorphic markers (e.g. Ali et al. 2014; Gladieux, Condon, et al. 2018; Spanner et al. 2021).

Host-specificity refers to the ability of certain pathogens to induce disease in particular host species or lineages within a host species (Borah et al. 2018; Li et al. 2020). The successful infection of a particular host species by a pathogen relies on the ability to overcome or avoid general and specific immune responses (Thines 2019). Molecular studies have provided mechanistic insights into the molecular interactions of plants and pathogens and have highlighted certain proteins on both host and pathogen sides, which are crucial determinants of resistance and susceptibility (reviewed in e.g. Dodds et al. 2006; Dodds and Thrall 2009). On the pathogen side, so-called effectors play determining roles in successful host invasion (Lo Presti et al. 2015; Toruño et al. 2016). Effectors are molecules secreted by pathogens to manipulate the host's physiology to the benefit of the pathogen. For example, the effector “EPIC1” is a protease inhibitor secreted by the oomycete pathogen Phytophthora infestans and related species infecting different plant species (Dong et al. 2014). EPIC1 targets immune-related and plant-specific proteases, thereby facilitating pathogen invasion in a host-specific manner. Another example is the effector Cmu1 produced by the biotrophic pathogen Ustilago maydis. This effector was shown to interfere with the salicylic acid biosynthetic pathways and thereby directly influence immune signaling (Djamei et al. 2011). Many described effectors are proteins, and the genes encoding these can be predicted according to specific sequence hallmarks. Effector genes are typically cysteine-rich (increasing their stability in the extracellular space) and contain a signal peptide for secretion into the plant apoplast or translocation into the host cytoplasm (Sperschneider and Dodds 2022). Given their fundamental role during pathogen infection, most effector-encoding genes are rapidly evolving in an ongoing evolutionary arms race with host targets (Möller and Stukenbrock 2017).

Host-specific lineages have been reported in many major fungal plant pathogens such as the rice blast fungus (Pyricularia oryzae, syn: Magnaporthe oryzae), for which cereal- and grass- specialized lineages were reported to be diverging despite the gene flow between them (Gladieux, Condon, et al. 2018). Other examples include the cultivar-driven divergence of the fungal wheat pathogen Zymoseptoria tritici (Hartmann et al. 2018) and closely related species of the oomycete pathogen Phytophthora infecting distinct host species (Dong et al. 2014). In most cases, however, the underlying traits that define host-specificity are not known and can only be predicted by comparative analyses of pathogens on different hosts.

The fungus Cercospora beticola is the causal agent of the disease Cercospora leaf spot (CLS) in cultivated beets (Beta vulgaris spp. vulgaris, which includes sugar beet, fodder beet, leaf beet, and garden or table beet, as well as wild sister species of these crops, including sea beet (B. vulgaris L. spp. Maritima (L.) Arcang.) (Tehseen et al. 2023). The pathogen undergoes multiple rounds of asexual reproduction within a single growing season (Rangel et al. 2020). Between sugar beet growing seasons, the fungus primarily survives in the form of hyphal structures called pseudostromata known to be resistant to desiccation (Rangel et al. 2020). These specialized structures, consisting of fungal tissue and remnants of host tissue, are situated within the leaf substomatal cavities of leaf remains and can remain viable for up to 2 years. While a known sexual form has not yet been identified, C. beticola is classified as a heterothallic ascomycete fungus. This classification arises from the identification of 2 oppositive mating-type idiomorphs (MAT1-1-1 and MAT1-2-1) in different isolates within populations (Groenewald et al. 2006). Previous studies found both mating types to be present at equal frequencies in certain regions of Europe and North America, providing indirect evidence for sexual reproduction (Bolton, Secor et al. 2012; Bolton et al. 2014; Vaghefi et al. 2017). Despite the wide host range of the pathogen, including on wild and domesticated plants, information about the evolutionary history of the pathogen is scarce, and it is not known if the wide host range of the fungus in fact represents different host-specialized lineages or a generally broad host range of individual isolates (Rangel et al. 2020).

Domesticated beets and one of the wild ancestors, the sea beet, have a long history associated with humans (Tehseen et al. 2023). The wild beet species B. vulgaris subsp. Maritima has been described as highly resistant to a variety of pests and pathogens and thereby as a potential source of new disease resistance traits for sugar beet (Tehseen et al. 2023). However, little is known about the diversity and biology of pathogens in B. vulgaris subsp. Maritima. Archaeological and historic evidence suggests that the wild ancestor of beets originates from the Mediterranean region (Biancardi et al. 2020). The first reports of beet cultivations came from the ancient Greeks and Romans, who were cultivating beets for their leaves. Beetroot, the part we commonly associate with beets today, was only developed during the late stages of the Roman Empire (Biancardi et al. 2020). By the 16th century, beet cultivation had spread to Northern Europe, and fodder beets were widely used for animal feed. In the 18th century, a German chemist, Andreas Marggraf, discovered that sugar could be extracted from beets, laying the foundation for sugar beet production (reviewed in Rangel et al. 2020). This marked a significant turning point, as beets became a vital source of sugar during the Napoleonic wars when sugar cane became scarce. Eventually, human migration to North America further dispersed cultivated beets during the last few centuries. Altogether, there has been a rapid diversification and a worldwide dispersal of beets associated with human history and crop use.

The global spread of plant species is often accompanied by the undesirable spread of associated pathogens (Stukenbrock et al. 2006; Sotiropoulos et al. 2022). Pathogen dispersal can lead to new adaptations allowing pathogens to establish in new environments, eventually in agricultural environments with distinct crop cultivars or management strategies, e.g. locally used fungicides (Croll and McDonald 2017). Agricultural environments may pose very strong selection pressures on pathogens, and evolutionary analyses of genetic data may reveal genetic loci that have been important for the establishment of the pathogen in certain environments. Evidence of such adaptation can sometimes be detected as “selective sweeps” in the genome. Typical signatures of selective sweeps are reduced genetic variation and increased linkage disequilibrium (LD) across a genomic region spanning the advantageous allele (Stephan 2016). A Click or tap here to enter text. number of methods are available to detect such signatures in population genomic data by integrating different signals in the data, for example, the local site frequency spectrum and LD patterns, and accounting for an underlying demography model (Pavlidis et al. 2013; Alachiotis and Pavlidis 2018). We previously demonstrated that selection for fungicide resistance has resulted in strong selective-sweep-signatures in North-American populations of C. beticola (Spanner et al 2021). We hypothesize that population genomic analyses may allow the further identification of functionally relevant traits that have been important for recent adaptation.

In the present study, we analyzed the genomes of 326 C. beticola isolates obtained from diverse beet fields and natural sea beet habitats from different geographical locations. Our objective was to infer the global population structure of C. beticola and identify the factors shaping this. In particular, we aimed to assess the genetic footprints of host specialization and local adaptation within the species. We utilized single nucleotide polymorphism (SNP) data to explore genetic variation in C. beticola. Our findings strongly support a host-driven divergence of pathogen populations. Furthermore, we scanned the genome alignment to detect signs of recent selective sweeps and observed a significant overlap between these signatures and potential virulence factors and fungicide resistance alleles. Our study reveals host-driven divergence of an important plant pathogen but also highlights the need for surveillance of pathogens on wild crop relatives.

Results

Generation of C. beticola Population Genomic Datasets

To explore the genetic structure of C. beticola populations and shed light to the emergence and recent history of this economically important sugar beet pathogen, we compiled a dataset encompassing whole-genome sequence data of 326 isolates from cultivated sugar beet across 4 continents (Africa, America, Central Asia, and Europe) (supplementary table S1, Supplementary Material online). The isolates from the UK were sequenced in Germany at the Max Planck Genome Center Cologne, while the rest of the isolates were sequenced in the USA at Novogene (Sacramento, California, USA) with Illumina sequencing technology. All sequencing reads were mapped to the reference genome of C. beticola (De Jonge et al. 2018) to call SNPs. The read coverage across the 326 genomes ranged from 11× to 32× with an average of 21×. We identified a total of 1,068,644 SNPs among the 326 isolates. A summary of the mapping statistics and SNP calling is summarized in supplementary table S1, Supplementary Material online.

Increased Genetic Diversity in Local Mediterranean Populations of C. beticola

The level of standing genetic variation in populations can give insight into the recent and ancestral demographic history. To explore genetic variation in C. beticola, we first conducted a comparative analysis of the genetic diversity among isolates from the 16 geographical locations, revealing significant differences in median values of the parameter π, a measure of overall nucleotide diversity within populations (P-value < 2.2e−16, df = 16, Kruskal–Wallis test). Subsequently, we performed Wilcoxon tests to assess significant differences between pairs of populations (Fig. 1a, supplementary table S3, Supplementary Material online). Our analysis showed that the 2 Mediterranean populations, Spain and Croatia, comprise significantly higher levels of nucleotide diversity (π_Spain = 0.0022, π_Croatia = 0.0020) than every other population. We consider that the higher diversity in the Mediterranean populations could reflect that the origin of the pathogen is in Southern Europe.

Also, the populations originating from wild sea beet plants collected at the sea shore in the UK at Southwold and Bawdsey showed higher genetic diversity (mean π_Swd = 0.0016, π_Bawd = 0.0014) compared with other C. beticola populations. Interestingly, the UK population isolated from sea beets in Orford showed the lowest genetic diversity of all C. beticola populations (mean π: 0.0008), suggesting a considerable variation in diversity among geographically close populations on the wild sea beet host. The Orford population furthermore showed a negative Tajima's D value (Fig. 1b) reflecting an excess of low-frequency alleles possibly resulting from a population expansion after a recent bottleneck. We speculate that the pathogen population at Orford has a different demographic history compared with the other pathogen populations on UK sea beet. The different levels of genetic variation moreover indicate a “metapopulation structure” with limited gene flow between C. beticola populations on sea beet in the UK. Here we define “metapopulation” as genetically differentiated subpopulations with very limited gene flow but which occur in close geographical proximity.

Overall, Tajimas D for each of the populations was close to 0, indicating a mutation-drift equilibrium in the populations. However, some differences in Tajimas D values are apparent. We compared Tajima's D values for the other geographical C. beticola populations using the Kruskal–Wallis test, which revealed significant differences among pairs of populations (P-value < 2.2e−16, df = 16, Kruskal–Wallis test). Subsequently, we performed a pairwise Wilcoxon test to assess significant differences between the populations (Fig. 1b, supplementary table S4, Supplementary Material online). In general, sea beet-infecting populations show higher levels of Tajima's D, which may reflect distinct demographic processes acting on these populations. The C. beticola populations in Spain, Denmark, and Orford showed significantly different Tajimas D values (D < 0) compared with the rest of the populations. We speculate that an excess of low-frequency alleles (reflected by the negative Tajimas D values) may reflect a recent bottleneck in Denmark and Orford, as the level of nucleotide variation in these populations likewise is low. On the other hand, the high nucleotide variation in Spain and the negative Tajima's D value may reflect a strong population expansion.

We next addressed the mode of reproduction in the pathogen populations, considering that C. beticola can reproduce either asexually or by mating between individuals with compatible mating types, Mat1-1 and Mat1-2. To this end, we determined the relative frequencies of mating type idiomorphs and applied a Chi-squared test to assess statistical significance of the deviation from a 1:1 ratio, expected in a population with random mating. This analysis revealed equal mating type ratios for all populations, i.e. consistent with sexual reproduction, except for the C. beticola population from wild beet in Orford (Mat1-1-1:Mat1-2-1 = 3.5, Chi-squared test, P = 0.001, df = 1) (Fig. 1c, supplementary table S4, Supplementary Material online). In the case of the Orford population, the Mat1-1 idiomorph is not present, suggesting a primarily asexually reproducing population. Overall, however, our mating type analysis provides evidence for a prominent role of sexual reproduction in C. beticola populations across the different geographical regions and on the different hosts represented with our sampling.

We also addressed recombination at the genomic level. We computed the extent of LD decay for each population. LD patterns can be informative about the extent of sexual reproduction. Recently founded populations, as well as, populations that primarily reproduce asexually, will show longer blocks of LD compared with recombining populations. In line with the analysis of mating types, we found significantly longer linkage blocks in the Orford, UK and Danish populations. The LD statistic r² was reduced to half its maximum value at 20 Kbp for the Danish population, while it never dropped to half for the Orford population (Fig. 1d, Table 1). Considering that the Danish population showed a balanced mating-type ratio, we propose that the long LD blocks derive from a recently founded, sexually reproducing population, while for the UK, Orford population, we attribute the long LD blocks to asexual reproduction.

The Host-driven Population Structure of C. beticola

We next characterized the population genetic structure of C. beticola based on a combination of methods utilizing the “independent SNP” dataset (see Materials and Methods). Firstly, we assessed the degree of clustering by implementing a detrended correspondence analysis (DCA) (Hill and Gauch 1980) (Fig. 2a and b). The DCA plot primarily reflected clustering according to the host species, and separated isolates obtained from farmlands and from the wild hosts. Interestingly, the isolates coming from sea beet and table beet clustered together. We further explored population structure using an ADMIXTURE analysis to infer the proportions of shared ancestries among the isolates (Alexander et al. 2009). We run the analysis for a number of clusters (K) ranging between 2 and 9 (supplementary fig. S1, Supplementary Material online). The number of clusters corresponding to the deepest subdivision was selected, using a cross-validation analysis. Models with a K value larger than 6 induced only a small decrease in the cross-validation error (Fig. 2b, supplementary table S5, Supplementary Material online). At K = 2, isolates from UK-sea beet and New York table beet shared the same ancestry, while the sugar and fodder beet-infecting isolates constituted the other cluster. At K = 3, the isolates were further separated within the sugar-fodder beet cluster, revealing 2 globally occurring sugar beet infecting clusters. Finally, at K = 6, we observed further sub-clustering within the UK sea beet-infecting isolates, mainly driven by the geographic origin of the isolates (at Bawdsey, Southwold, and Orford). In addition, at K = 6, the North Dakota sugar beet-infecting population was divided into 2 genetic clusters, although obtained from the same geographical location. We note that shared ancestry, in general, is extensive among the C. beticola isolates investigated here. In line with the analysis of genetic diversity, the Croatian population encompasses most of ADMIXTURE's hypothetical ancestral clusters, possibly reflecting that the population is older, and results from a longer history of admixture of different populations in the region.

Host Species is a Significant Factor Separating C. beticola Populations

Host specialization and geographic isolation are 2 factors that often play a key role in reproductive isolation in filamentous plant pathogens (Möller and Stukenbrock 2017). To illuminate the role of geographic separation and host specialization as the driving factors of C. beticola population divergence, we first performed 2 AMOVAs (supplementary fig. S2, Supplementary Material online, Table 2). To assess the fraction of variation associated with the geographic distance, we used all the available isolates. In this analysis, the major component of genetic variation, 77%, is explained by variation within individual populations. An additional 19% of the variation is explained by genetic variation within countries. Only 9% of the variation is explained by variation between countries. Furthermore, we found did not find evidence that geography at the continental scale explains variation between C. beticola populations. Overall, the AMOVA suggested that genetic variation within populations explains most of the variation in C. beticola, while variation between countries and continents is nonsignificant (P-value: 0.17 and 0.82, respectively). This observation indicates a recent dispersal of C. beticola and/or extensive gene flow across long distances. We further generated a neighbor network, which further supports a lack of geographical population structure and a recent differentiation of the populations (Fig. 3). The neighbor network moreover suggests that genetic variation has been homogenized by gene flow among isolates from different countries.

We next tested the effect of the host on the population structure of C. beticola. For this analysis, we focused on the UK, the only country from which we had populations obtained from different hosts within the same region. Although, the main contribution of the variance again is explained by variation within populations, the effect of the host species (sugar vs. sea beet) was significant (AMOVA, P-value: 0.03, df = 1). As much as 21% of the genetic variance was associated with the host, suggesting a strong factor of host genetics in driving genetic differentiation of the UK C. beticola populations. We note that other factors, such as local adaptation and genetic drift, may further contribute to the population structure of the pathogen. Altogether, we conclude, however, that for the geographically coexisting C. beticola populations on wild and cultivated beet in the UK, there is a significant effect of host type on the partition of genetic variance.

Phylogenomic Analysis Suggests an Ancient Split of Wild Beet-infecting C. beticola Populations and Recurrent Introductions to the Agro-ecosystem

Domesticated beet was dispersed across Europe during the medieval era and later to the Americas with European migrants (Biancardi et al. 2020). To test if early dispersal of C. beticola coincided with the dissemination of beet cultivation, we assessed the evolutionary relationships of the C. beticola populations. For this purpose, we used a maximum likelihood phylogenetic analysis, IQ-TREE, implementing a polymorphism-aware model (PoMo) (Minh et al. 2020) (Fig. 4). Using the sister species Cercospora flagellaris as the root for the C. beticola populations, the analysis revealed 2 primary population clades: a lineage comprising the Mediterranean, sugar, and fodder beet-infecting populations, and 1 lineage comprising the UK sea beet-infecting isolates and the NY Table beet-infecting isolates. An intriguing conclusion from the phylogenetic analysis is that C. beticola was thereby introduced at least 2 times to agricultural ecosystems. We base this conclusion on the fact that the 2 main clusters each comprise a distinct sea beet population, 1 in Southern Europe (Croatia), which may have given rise to the sugar beet-infecting C. beticola pathogen, and 1 in Northern Europe (the UK sea beet population), which may have given rise to the table-beet-infecting C. beticola population.

Ongoing Host-driven Divergence of C. beticola in the UK

Despite the fact that the UK isolates obtained from farmlands and wild plants were collected from a distance of roughly 150 km apart, our ADMIXTURE analysis suggests that they have different ancestries (Fig. 2). Furthermore, there was no evidence of admixture between these populations. Indeed, a Principal Component Analysis (PCA) including only the UK populations showed a clear separation between the sugar and sea beet-infecting populations (Fig. 5a, supplementary fig. S3, Supplementary Material online). As the AMOVA described above indicates, some genetic variation is explained by the host from which the isolates were obtained. Therefore, we hypothesized that divergent host specialization may act as a barrier to gene flow. To pinpoint the genomic components associated with host specialization in the UK, we performed genome-wide scans for highly divergent regions between the 3 sea and 1 sugar beet infecting populations sampled in the UK. We computed the population divergence statistic F_ST and D_XY for nonoverlapping windows of 5 Kbp along the genome. We considered genomic windows showing F_ST and D_XY values in the 0.05%ile to be outlier genomic regions representing loci of increased divergence between populations. To further control for the effects of varying genetic diversity that can influence the F_ST results, we focused on the regions that have been identified as significant by both statistics. A map of the highly diverging regions between each of the 3 sea beet and the sugar beet populations is shown in Fig. 5b.

Due to their determining role in host-pathogen interactions, genes encoding effector proteins often experience strong selective pressure (Möller and Stukenbrock 2017). Effector proteins play a crucial role in plant-pathogen interactions, as pathogen lineages with the right effector repertoire can suppress the defenses of a specific host and successfully complete their life cycle. Therefore, we conducted a permutation test to determine whether effector genes showed a significant overrepresentation in the highly diverging regions. To this end, we evaluated the proportion of predicted effector genes within the highly diverging regions in comparison to the remainder of the genome. Our results clearly indicated that effector genes are enriched in the highly diverging regions when compared with randomly selected genomic regions, with 62 predicted effector genes in these regions (Fig. 5c). The over-representation of effector genes in highly diverged regions supports our hypothesis that this group of genes is under strong selection, likely due to their prominent function in virulence. We speculate that different repertoires of immune receptors and resistance genes in beet cultivars and species drive adaptive changes in effectors of C. beticola.

The Recently Described Avirulence Gene, AvrCR4, has Been Present at Low Frequencies in C. beticola Populations

Sugar beet lines harboring the major resistance gene BvCR4 are resistant to C. beticola (Chen et al. 2024). This resistance gene was recently identified and introduced into cultivated sugar beet. The BvCR4 resistance gene has, however, rapidly lost its efficacy as virulent C. beticola have emerged in the last few years. In a recent publication, the avirulence gene was identified and described as an effector-like small, secreted protein (Chen et al. 2024). Deletion of the avirulence gene AvrCR4 is associated with high levels of virulence on BvCR4 beet lines. With the available sequence information of AvrCR4, we scanned our geographical population samples for individuals that carried the AvrCR4-deletion. In line with the findings from Chen et al. 2024, we find that the deletion is very rare, and observed in only 2 out of the 289 isolates we analyzed in this study. The 2 sugar beet isolates that carry the deletion originate from the UK and Sweden. Intriguingly, our samples were collected between 2011 and 2017, before the release of the widely used BvCR4 sugar beet variety into the market in 2020. Therefore, the recorded frequency of the deletion reflects the period before the introduction and shows that the deletion mutation was already present in C. beticola populations before the introduction of the BvCR4 resistance gene.

High Frequencies of Fungicide Resistance-associated Alleles in C. beticola Populations

In farmlands, there is often an extensive use of pesticides, such as fungicides applied to manage fungal disease (Stukenbrock and Mcdonald 2008). Demethylation inhibitor compounds (DMIs) are prominently used to control C. beticola populations in sugar beet fields worldwide. In a previous study, 4 SNP alleles were reported to be associated with DMI resistance (Spanner et al. 2021). To examine the difference in the abundance of the identified DMI resistance alleles, we recorded their frequency in the 4 UK populations, including isolates from wild and cultivated beet (Fig. 5d). The mutations associated with fungicide resistance were located at positions 2,402,041 on chromosome 1, 849,506 on chromosome 4, and 1,358,331 and 1,451,478 on chromosomes 9. All the mutations are located in hypothetical proteins except the latter, which is a synonymous mutation within the coding region of the gene encoding for eburicol 14-alpha demethylase CB0940_11379, otherwise known as the gene encoding the DMI fungicide target CbCYP51 (Bolton, Birla, et al. 2012). Three out of 4 mutations were absent in the sea beet-infecting populations, while the mutation on chromosome 4 was found in all populations at varying frequencies (15% to 33%). In contrast, the mutations on chromosome 4 and 9 were found in the UK sugar beet infecting populations in high frequencies (in 75%, 53%, and 53% of all isolates for chromosome 4, chromosome 9: 1,358,331, and 1,451,478, respectively). None of the UK populations carried the mutation on chromosome 1. Both mutations on chromosome 9 are in close vicinity of genomic regions that we found to be significantly diverged between the UK farmland and wild populations (Fig. 5e). Altogether, this analysis shows a strong selection for DMI fungicide resistance in UK farmlands, while these mutations are rare or absent in populations infecting wild beet plants. We further recorded the frequencies of DMI fungicide resistance in C. beticola populations worldwide. Fungicide resistance-associated variants were present in all populations around the world at different frequencies (supplementary figs. S4 to S14, Supplementary Material online), and we conclude that genetic variants associated with DMI resistance are present in both the farmland and wild populations of C. beticola and may be readily exchanged between isolates.

Discussion

Lineages of C. beticola Have Recently Specialized to a New Crop Species in Agro-ecosystems

Understanding the factors driving the evolutionary trajectories of crop pathogens is essential to designing surveillance strategies and preventing future epidemics. Here, we analyzed a population genomic dataset of the globally present pathogen of sugar beet, C. beticola, to investigate patterns of genomic variation and signatures of divergent host specialization. To this end, we used an extensive global sampling of the pathogen from wild and domesticated hosts to assess the population structure and the extent of host specialization on wild and domesticated plants. Our analyses were designed to test the hypothesis that agricultural practices and the domestication of beets, involving substantial genetic changes and phenological differences between wild and cultivated plants, have shaped the recent evolution of C. beticola lineages. Our detailed population genomic analyses provide evidence that C. beticola is subdivided into multiple lineages associated with wild and domesticated hosts. Our inference of ancestry contribution revealed 6 main ancestral clusters. Intriguingly, isolates from the wild showed different ancestries than those from farmlands. This separation was even more surprising for isolates collected at a relatively small distance (within 200 km) in UK. The UK isolates from farmlands and wild sea beet showed different ancestries, and there was no evidence of extensive gene flow between the populations. This indicates some extent of host specialization among closely related lineages of C. beticola. Host specialized lineages were also reported for major plant pathogens like the rice blast fungus (Pyricularia oryzae, syn: Magnaporthe oryzae) (Gladieux, Condon, et al. 2018) and the wheat pathogen Zymoseptoria tritici (Hartmann et al. 2018). In Hartmann et al. (2018), the authors used statistical methods to identify genes under selection in different host-specific lineages of the pathogens. Signatures of positive selection were identified in known effector genes, indicating adaptation to distinct wheat cultivars. Moreover, genes associated with fungicide resistance were found to co-localize with selective sweep regions in the wheat pathogenic fungus Z. tritici (Hartmann et al. 2018). Highly divergent genomic regions between C. beticola populations from the wild and from farmlands were also found here to be enriched with predicted effector genes.

Further evidence of host-specific divergence was obtained from an AMOVA. AMOVA results showed that the effect of geographic distance is not significant on the partition of genetic variance of the species. That is in line with previous works that assessed the population structure of C. beticola isolates from farmlands and concluded the presence of a shallow population differentiation (Vaghefi, Kikkert, et al. 2017; Vaghefi, Nelson, et al. 2017; Knight et al. 2019). Our analysis included isolates of C. beticola from the wild ancestor of beet, the sea beet. This allowed us to assess the significance of the host in the species population structure and estimate the proportion of genetic variance associated with it. The analysis showed that the host significantly explained 21.42% of the variation between C. beticola populations and supported the hypothesis that the different hosts drive genetic differentiation in pathogen populations.

Sugar Beet Infecting Populations of C. beticola Have Likely Emerged in the Mediterranean Region with Subsequent Introgressions From Wild Beet Species

A primary indicator of the center of origin of a species is the high genetic diversity compared with other locations where the species was later introduced. In the case of C. beticola, we observed higher diversity primarily in the Mediterranean region and secondarily in the wild sea beet population collected in the UK. Some studies have highlighted the effect of host domestication in the emergence of its associated pathogens, including the wheat-infecting pathogen Z. tritici responsible for the Septoria tritici blotch (Stukenbrock et al. 2007), the rice blast fungus Magnaporthe oryzae (Couch et al. 2005), the corn smut fungus Ustilago maydis (Schweizer et al. 2021), and the barley net blotch fungus, Pyrenophora teres (Taliadoros et al. 2024). The Mediterranean Sea was previously proposed as the centers of diversity and origin of cultivated beets (Biancardi et al. 2020). An elevated diversity in that region indicates that the populations are older and may represent a center of origin. We speculate that the emergence of C. beticola could have coincided with the early domestication of the beet host in Southern Europe.

In our study, we also assessed the evolutionary relationships between the populations using polymorphism-aware phylogenetic models. We found a common origin of the sugar beet and fodder beet infecting populations that formed a clade with the Croatian wild population, which was early diverging within the clade. This suggests a Mediterranean origin of the sugar/fodder beet populations. In the second clade, the UK wild and NY table beet populations were clustered, which suggests that the UK origin of the NY population has a different origin than that of the other pathogen populations included in this study. The different origins of the sugar beet-infecting and table beet-infecting populations lead to an intriguing conclusion that C. beticola was introduced to the agro-ecosystem on multiple occasions. Based on the widespread occurrence of ancestral variation, we speculate that the migration history of C. beticola reflects a recent range expansion of different ancestral genotypes giving rise to present-day admixed populations.

Our extended C. beticola dataset allowed us to conduct a detailed exploration of the demographic history of the species and assess the impact of host domestication, which should be a future step to further understand the emergence and spread of the CLS disease.

Sexual Reproduction Contributes to Evolution of C. beticola Populations

We used SNP data to calculate genetic diversity within C. beticola populations. Interestingly, C. beticola exhibited a higher level of nucleotide diversity (π) in comparison to well-studied crop pathogens like the wheat pathogen Z. tritici (Hartmann et al. 2018), the wheat powdery mildew pathogen B. graminis f.sp. tritici (Sotiropoulos et al. 2022), and the rice blast fungus M. oryzae (Zhong et al. 2018). Genetic variation plays a crucial role in the ability of a pathogen to rapidly adapt; therefore, the elevated genetic diversity in C. beticola may be a significant factor contributing to the pathogen's successful survival and dispersal. Disparities between species may be attributed to variations in the frequency of sexual recombination, historical demographic events, such as bottlenecks, and gene flow. Our analyses of LD and the mating-type frequency distribution support the hypothesis that C. beticola frequently undergoes sexual recombination. This is in line with previous reports for North American and European populations (Vaghefi, Nelson, et al. 2017; Knight et al. 2019). An exception is observed in the UK population from Orford, where clonality appears to be more prevalent. Alternatively, the Orford population represents a very young founder population in which genetic diversity is still to be built up. The population genetic structure of the UK population reveals the presence of 4 distinct genetic clusters, 3 obtained from natural environments and 1 from sugar beet fields. Despite these isolates coexisting in a relatively small region, they belong to distinct clusters with limited signs of gene flow between them. This observation suggests potential host specialization, reflecting a meta-population structure and constraints on gene flow, particularly between isolates adapted to farmlands and to wild plants. A recent study of Z. tritici in Iran likewise found evidence for a metapopulation structure where populations of the pathogen infecting different hosts were genetically isolated (Fagundes et al. 2024).

Host Species and Agricultural Practices Drive Adaptive Evolution in C. beticola

To identify evidence of local adaptation, we utilized a variety of methods to identify genomic regions that are highly diverging between the 1 farmland and the 3 wild populations collected in the UK. In total, we identified 214 highly diverging regions between these populations. Diverging genomic regions between farmland and wild C. beticola populations could harbor important host-specificity loci, and we indeed find that these were enriched with predicted effector genes, further highlighting the importance of the host as a main driver of rapid evolution in this pathogen.

A recent study reported genomic loci associated with DMI resistance (Spanner et al. 2021). We recorded the abundance and frequency of DMI resistance alleles in the 4 UK populations. We found 3 of 4 variants associated with fungicide resistance present in high frequencies in the UK farmland population. In contrast, we found only 1 of the 4 variants in low to moderate frequencies present in the wild populations. Furthermore, 2 of the DMI resistance loci are located in the highly diverging genomic regions, possibly reflecting a rapid adaptation of the pathogen to the agro-ecosystem. DMI compounds are used to manage C. beticola populations in all sugar beet growing regions, including the UK (Liebe et al. 2023). Furthermore, we found fungicide-resistant associated alleles in high frequencies in various populations around the world.

Our data does not allow us to identify the origin and source of genetic variants linked to DMI resistance. However, we consider 2 likely scenarios: (i) The mutations may have already been present in the wild as part of natural genetic variation and subsequently increased in frequency in farmland populations in response to the selective pressures imposed by field treatments. (ii) Alternatively, these could be novel mutations that evolved within farmland populations of C. beticola and subsequently introgressed into populations of the pathogen on wild hosts. Under both scenarios, C. beticola populations on wild beet species may act as a reservoir of virulence or fungicide resistance-related alleles. Further research is needed to unravel the emergence and spread of DMI resistance in C. beticola and the relevance of genetic exchange between populations on wild and cultivated hosts. Our study underlines the importance of studying pathogen populations on wild plants in order to understand the rapid evolution of crop pathogens. Close relatives of crop pathogens that infect wild plant, may be sources of genetic variants that confer new virulence specificities or fungicide resistances, as demonstrated for the C. beticola pathogen.

Conclusion

Using a population genomics framework, we show that C. beticola is subdivided into multiple lineages associated with cultivated and wild hosts. We further present evidence of limited genetic exchanges between them. Our genetic diversity and phylogenetic analyses provide evidence of a parallel demographic history between the host and the pathogen and multiple introductions of the pathogen to the agro-ecosystem. This adds to the growing body of work highlighting the important role of human activities, such as domestication, migration, and trade in the emergence and spread of plant pathogens. Furthermore, our results enhance our understanding of the evolutionary forces driving the diversification of C. beticola and underscore the importance of research on wild plant pathogens as well as the utility of genomic data for epidemiological surveillance of fungal plant pathogens.

Materials and Methods

Genome Data

Out of 326 C. beticola whole genomes were sequenced using Illumina technology. The sequenced isolates were sampled from cultivated and wild hosts on 4 different continents (Table 1 and supplementary S1, Supplementary Material online). Eighty isolates were collected from North Dakota State University experimental sugar beet fields in Fargo and Foxhome, North Dakota, USA. Twenty-four isolates were collected from 2 table beet fields in New York, USA. Seven isolates were collected from a fodder beet field in Chile. Twenty-two isolates were collected from a sugar beet field, in the UK. Furthermore, 12 isolates were collected from a sugar beet field in Sweden. Isolates were also collected from single sugar beet fields in 6 continental European countries. Specifically, 12 isolates were collected from a field in Denmark, 18 from Germany, 11 from France, and 19 from Spain. Additionally, 21 isolates were collected from 5 fields in Italy. Moreover, we included 18 isolates from Turkey. To capture the natural variation of the species in continental and Northern Europe, we included 23 isolates from 2 natural habitats of sea beets in Croatia and 43 isolates from 3 natural sea beet locations in the UK (Southwold, Bawdsey, and Orford) (Table 1).

Read Mapping and Variant Calling

To extract high-quality SNPs, we processed and aligned Illumina sequencing reads to a reference genome. Initially, the sequencing adapters and low-quality nucleotide bases were filtered using the software Trimomatic version 0.38 (Bolger et al. 2014), applying a PHRED quality score of 33 and a minimum read length of 30 base pairs (bp). To control for potential overrepresentation artifacts, overlapping reads were merged with PEAR version 0.9.11 (Zhang et al. 2014).

Subsequent alignment of individual reads was performed using the Cb09-40 C. beticola reference genome (De Jonge et al. 2018). To this end, we used the Burrows-Wheeler Aligner version 0.7.17 (Li and Durbin 2010) and Stampy v. 1.0.20 (Lunter and Goodson 2011) (see supplementary table S2, Supplementary Material online for details). SNP calling was performed using the GATK HaplotypeCaller version 4.2.18 (McKenna et al. 2010), generating a variant call format (VCF) file containing raw SNP calls. Next, to remove low-quality SNPs, we followed a “hard-filtering” approach in accordance with GATK's best practice guidelines. Specifically, the following filtering criteria were applied: (i) the ratio of SNP call quality to sequencing depth had to be higher than 2, (ii) a minimum sequencing depth of 8 was required, (iii) SNPs had to have a mapping quality score >40, (iv) the allele-specific rank sum test for mapping qualities of reference (REF) versus alternate alleles had to be >−12.5, (v) the allele-specific rank sum test had to be higher than −8, and (6) each genome required an average read coverage of at least 2. The filtering process was performed using GATK VariantFiltration version 4.0.11 (McKenna et al. 2010).

After the hard-filtering process, a final dataset of 1,129,645 variants was obtained, 893,266 of which were biallelic SNPs. This dataset is referred to as the “full high-quality dataset.”

Many clustering analyses assume that the genetic markers are independent. To meet this assumption, the “full high-quality dataset” was further filtered based on patterns of LD decay, ensuring a minimum distance of 2 kilobase pairs (Kbp) between SNPs. This threshold was determined by calculating the decay of LD (r²/2) averaged across populations (Table 1). Following LD-based filtering, a final dataset comprising 44,229 biallelic SNPs was retained. This dataset is referred to as the “independent SNP dataset.”

Genetic Diversity, Neutrality Tests, and LD

We used the “full high-quality dataset” to compute and compare genetic variation among populations. The software pixy (Korunes and Samuk 2021) was applied to obtain the genetic diversity estimator π and the population divergence statistics F_ST and D_XY in nonoverlapping windows of 10 Kbp. The tool PopLDdecay (Zhang et al. 2019) with default parameters calculates LD decay. Tajima's D was estimated with VCFtools v. 0.1.17 (Danecek et al. 2011). Significant differences in the level of genetic diversity and Tajima's D values were assessed with the nonparametric Kruskal–Wallis test and the pairwise Wilcoxon with Bonferroni correction. The selection of these nonparametric tests, rather than a parametric test (e.q. ANOVA), was made because not all group data followed a normal distribution and the groups did not have equal variances, as parametric tests assume.

Population Genetic Structure

We inferred population genetic structure using 3 different approaches: we conducted a DCA (Hill and Gauch 1980), we applied an analysis of admixture using ADMIXTURE version 1.3 (Alexander et al. 2009), and we performed a Neigbor-Net analysis using SplitsTree v 6 (Huson and Bryant 2024). These analyses were based on the “independent SNP dataset.” DCA analysis was performed using the R package vegan v. 2.6.2 and DCA visualization with R package ggplot2 (Wickham 2016). We further characterized population structure using the software ADMIXTURE. We first ran ADMIXTURE assuming K-values between 1 and 10. The most fit K value was assessed based on the cross-validation error (Alexander et al. 2009) (supplementary table S5, Supplementary Material online). Finally, to address the effect of geographic distance and individual host on the population structure, we employed an analysis of molecular variance (AMOVA) (Excoffier et al. 1992) and performed this on the same dataset. To this end, the isolates were grouped geographically according to the country and the continent of their origin, and according to the host they were isolated from. The R package poppr (Kamvar et al. 2015) was used to perform the AMOVA. Since our dataset encompass isolates from different hosts collected in the UK, we could specifically assess the significance of geographical origin and host specialization in 2 independent analyses. To assess the significance of geography, we used data from all available C. beticola isolates, while for the analyses of the host specialization; we focused on the UK C. beticola populations from sea beet and sugar beet.

Mating Types

MAT1-1 and MAT1-2 mating-type sequences were obtained from GenBank (accession numbers DQ264736.1 and JN863091.1, respectively) (Groenewald et al. 2006; Bolton, Secor et al. 2012). SPAdes (Bankevich et al. 2012) with default parameters was run to create complete de novo assemblies of C. beticola genomes. Subsequently, BLASTN with default parameters (Camacho et al. 2009) was used to assess the mating type of each isolate.

Evolutionary Relationships Using PoMo

We employed the PoMo, as implemented in the IQ-TREE software (Minh et al. 2020) to reconstruct relationships among C. beticola populations. To convert our FASTA file to the count file format, we utilized the FastaVCFtoCount.py script provided with the PoMo software. We selected a substitution model based on the Bayesian Information Criterion. To this end, the ModelFinder algorithm implemented in the IQtree version 2.0.3 was used. To obtain branch support values, we perform 1000 bootstrap replicates, employing IQ-TREE bootstrap approximation option. The evolutionary relationship of the C. beticola populations was reconstructed with the sister species C. flagellaris as an outgroup.

Genome Scans for Population Divergence

The software pixy was used to estimate the population divergence statistics F_ST and D_XY in nonoverlapping windows of 10 Kbp. Genes encoding putative secreted proteins were identified as proteins exhibiting a signal peptide using the program SignalP (Almagro Armenteros et al. 2019). SignalP-5 predicts the presence of a signal peptide and the location of its cleavage sites. TMHMM v. 2.0 was used for the prediction of transmembrane proteins (Krogh et al. 2001), which were excluded for the further analyses. The remaining genes were used to predict putative apoplastic and cytoplastic effectors with the deep learning tool EffectorP v.3 (Sperschneider and Dodds 2022).

Next, the prevalence of predicted effector genes in the diverging regions compared with other genomic regions was examined using a permutation test (custom script available on GitHub: https://github.com/Jimi92/). In brief, the abundance of predicted effectors was counted in regions of the same size and number equal to the size distribution of the diverging regions. As many as 10,000 replicate runs of random resampled regions were performed. Genome-wide maps of the diverging regions were created using Circos v. 0.69-9 (Krzywinski et al. 2009).

Supplementary Material

Supplementary material is available at Genome Biology and Evolution online.

Acknowledgments

The authors are grateful to all members of the Environmental Genomics group. Furthermore, the authors acknowledge the British Beet Research Organization for access to trial sites.

Funding

D.T.'s work was supported by the DFG Research Training Group RTG 2501. This work was furthermore supported by the U.S. Department of Agriculture, Agricultural Research Service through In-House Appropriated Project 3060-21000-045-000-D.

Data Availability

Raw reads for the C. beticola isolates sequenced for this work are available in the NCBI short- read archive under BioProject PRJNA1174990. Moreover, information files, custom scripts and workflows are available at the Github repository: https://github.com/Jimi92/Cercospora_beticola_pop_gen.

Literature Cited