The analyzed population comprised 8,316 winter wheat (T. aestivum L.) PGR maintained at the Federal ex situ gene bank for Agricultural and Horticultural Crop Species, Germany, hosted at the Leibniz Institute of Plant Genetics and Crop Plant Research (IPK), as well as 171 European winter wheat varieties representing the currently cultivated varieties in Central Europe (in the following denoted as the Elite panel). For each PGR, seed samples were provided by the gene bank and multiplied in a first step using two-row plots (for details, see Schulthess et al., 2021). At this time, one and in rare cases two representative ears per plot were isolated during flowering for controlled self-fertilization and harvested separately from the rest of the plot at maturity. Seeds produced from each selfed single ear were then multiplied in a single-row plot during the following crop season. Seeds harvested from the single-row plots corresponded to the genetically defined propagation material of PGR used for DNA extraction and phenotyping. For the Elite panel, seeds for these purposes were collected from the seed market or directly provided by the seed industry.

The phenotypic data was gathered using the Macrobot facility, a robotic platform for performing high-throughput semi-automatic detached leaf assays developed at the IPK Gatersleben (Lück et al., 2020a,b). Seedlings for the Macrobot assay were grown on trays comprising 6 × 4 slots in a greenhouse with standardized conditions at the IPK Gatersleben. For each genotype, 10 seeds were sown in one slot. Then, 15 days after sowing, the second leaves were harvested from seedlings. Up to eight leaves of different seedlings of the same genotype were collected. Leaf segments with a uniform length of 2 cm obtained from the middle part of the harvested leaves were mounted onto 4-well microtiter (MTP) agar-plates (1% water agar supplemented with 20 mg L^-1 benzimidazole as a leaf senescence inhibitor), following the sowing pattern of the growth trays (see Figure 1 for a graphical illustration). The plates were inoculated in an inoculation tower with the highly virulent Blumeria graminis f. sp. tritici isolate FAL 92315 (see Supplementary Table 5). Spores from heavily PM infected plants were applied from the top of the inoculation tower using compressed air while the rotation of the platform with the open plates assured an even distribution of the spores. The inoculation density was controlled on microscope slides added between the plates and counted under a microscope after each portion of spores until the optimal spore density was reached (5–10 spores per mm²). The size of the rotatory platform allowed 12 MTPs with leaves harvested from two growth trays to be simultaneously inoculated, thus forming an experimental group. After inoculation, the plates were covered with the lids, transferred to an incubation chamber, and incubated for six days under controlled conditions (20°C, 60% RH constant, 16 h light μE m^-2 s^-1). Following incubation, multimodal images of the plates were automatically acquired and analyzed as described in Lück et al. (2020a) and Lück et al. (2020b). The disease severity was calculated as a percentage of the infected leaf area. In total, 8,487 genotypes were tested in 446 experiments. Each genotype was tested in two different experiments and up to eight individuals per genotype in each experiment. The susceptible cultivar KANZLER was used as running control by placing it in four out of the 24 slots of each growth tray.

Detached leaf assays are artificial systems that are prone to certain inconsistencies due to the complexity of pathogen-host interactions and environmental influences, which can be reduced but never completely eliminated. Manual assessment of data quality was not feasible given the volume of the data, so an automatic, standardized quality pipeline was developed. To ensure high data quality, unreliable data at technical replicate, experimental, and biological replicate levels were identified and excluded from further analyses. The following data curation steps were performed: The first stage of outlier correction was based on the technical replications. In some cases, border-effects or leaves in a generally unhealthy state resulted in extreme values within the 3–8 technical replications. These were excluded by using 1.5 times the interquartile distance as the threshold. In a second step, the data were analyzed at the experiment level. Based on the susceptible control KANZLER, experiments with generally very low infection were excluded. Outliers were defined by using the 1.5 interquartile distance again as a threshold of the control genotype alone. If the mean or maximal value of the control of an experiment was below this threshold, the experiment was excluded from further investigation. The final quality check was based on the variance between the biological replicates, that is replicates of the same genotype in two different experiments. We implemented the outlier detection based on residuals as described in Anscombe and Tukey (1963) considering a nominal α level of 0.05. For this, we fitted the same model as used in Equation (1) and extracted the residuals. Datapoints exceeding the significant threshold were flagged as significant outliers and were excluded from further investigation. All computational methods were performed within R environment (R Core Team 2018 version 4.0.2. using R-Studio version 1.3.1056).

Variance components of the phenotypic data were estimated using a linear mixed model approach (Henderson, 1975).

The following mixed model was fitted to the percentage of infected leaf area data (y):

where the common mean μ was treated as a fixed factor, whereas genotype, experiment, the tray nested within an experiment, and error effect were assumed as random factors. Best Linear Unbiased Estimations (BLUEs) of each genotype were computed using the same model as in Equation (1) but assuming the genotype factor as fixed. For variance components estimation and BLUEs computation, linear mixed models were solved using the ASReml-R package Version 4 (Bultler et al., 2009). The heritability was estimated as in the following equation:

where σG2 is the genotypic variance, σe2 is the residual variance while R represents the average number of replications per genotype. Heritability estimations from 500 different random samples containing 80% of the total number of genotypes were used to estimate the standard deviation of the heritability.

The study is based on genomic data presented in detail elsewhere (Schulthess et al., 2021). Briefly, genotyping was carried out at IPK facilities according to Poland et al. (2012) with some modifications. DNA was extracted from seedlings using established protocols by KWS LOCHOW GmbH. The DNA from each genotype was digested using PstI and MspI restriction enzymes (New England Biolabs) following Wendler et al. (2014) and ligated with adapters containing sample-specific barcode sequences. Later, processed barcoded DNA samples were pooled into groups of 540 genotypes in an equimolar amount to form a GBS library. Single-end (100 bp) sequencing was performed on Illumina HiSeq 2500. On average, 2.5 million reads were generated for each genotype, and these sequence reads were complemented with sequence information for 171 additional wheat varieties. Raw reads were trimmed in order to remove adapters and low-quality base calls. High-quality reads were then mapped against the reference genome sequence of Chinese Spring (IWGSC RefSeq v1.0), and SNP calling was carried out resulting in 1,646,929 SNP. The SNPs obtained were further filtered to keep those with homozygous genotype counts for each allele greater than 10, missing value rates smaller than 10%, and heterozygosity levels smaller than 1%. After filtering, information for 29,183 SNPs across 8,070 genotypes was available. To avoid introducing non-existent genetic variation into each pool, missing values were imputed separately for PGR and the Elite panel, respectively, using Random Forest implemented in the missForest R-package (Stekhoven and Bühlmann, 2012). Among genotypes with GBS profiles, 7,337 PGR-Isolates and 154 elite varieties had also curated Macrobot data available.

Genome-wide association mapping was performed using the R-Package rr-BLUP (Endelman, 2011), which uses a mixed model according to Yu et al. (2005):

where y is a vector of the phenotypic data (BLUEs of infected leaf area); β is the fixed population mean; α is the fixed effect of the tested SNP; u is a random vector of polygene background effects and e is a vector of residual effects. X, S, and Z are incidence matrices relating y to β, α, and u. To account for population stratification and relatedness, a kinship matrix based on two times one minus Rogers' distance (RD) (Rogers, 1972) was estimated. Linear mixed models for GWAS were solved by using the “efficient mixed-model association” (EMMA)-algorithm (Kang et al., 2008) as implemented in the rrBLUP R-package (Endelman, 2011). Two significance thresholds were used for multiple-testing correction, both using a nominal level of test significance α = 0.05. Both thresholds were obtained using Bonferroni correction (α/n) but with two different n values. For the more restrictive threshold, n was defined as the total number of markers (29,183), while n of the less conservative threshold was calculated as the effective number of independent markers, M_eff (5,996) following Gao et al. (2008).

For each significantly associated SNP, we compared the subpopulations defined by the two different alleles. The allele that occurred more frequently in the more resistant group was determined to be the beneficial allele. The significant SNP-trait associations were further analyzed using the Reference Genome Chinese Spring (IWGSC RefSeq v1.0). In a first step, we combined adjacent significant loci that were not interrupted by non-significant markers, aggregating marker-trait associations (MTA) to associated genomic regions.We also calculated the linkage disequilibrium (LD) for the SNPs within a defined region as the squared correlation coefficient (r²). We then collected the information of the predicted high-confidence genes and their predicted functional domains within these regions using the EnsemblePlants platform (Howe et al., 2019). The lines carrying the beneficial allele were defined as possible donor genotypes. We estimated the (r²) for the loci by fitting them individually into a simple linear model.

To test whether certain domains were enriched around the significant associations, a domain enrichment analysis was performed. Using the Pearson's Chi-squared test, this analysis compares the relative abundance of the domains predicted in defined regions and their relative abundance in the whole genome. A significant region is defined by the distance between the next non-significant marker downstream and upstream of every significant marker-trait association. For this analysis, the reference genome annotation IWGSC v1.0 HC 20170706 was used. This version of the reference genome annotation contains 105,589 high confidence gene predictions. For all those genes, we predicted domains using PANTHER (Thomas et al., 2003), Gene3D (Lees et al., 2012), Pfam (Mistry et al., 2021), and PROSITE (Sigrist et al., 2012). First, for each domain type, the relative abundance was calculated. Then the relative abundance of each domain was calculated only using genes located in one of the significant regions.

Phenotypic data were compared with data from field trials of the German Federal Plant Variety Office (Bundessortenamt, BSA). In these trials, the severity of the PM infection was rated by visual evaluation in nine categories from 1 (no infection) to 9 (severely infected). The data set contained phenotypic data from 2002 till 2019 for a total of 365 wheat varieties. Among these lines, 109 were also tested with the Macrobot. The Best Linear Unbiased Estimations (BLUE) were obtained by fitting the following linear mixed model:

where, y is the grade of infection, μ and genotype are the fixed general mean and fixed genotype effects, respectively, while year is the random effect of the year.

Powdery mildew resistance screening of the winter wheat collection using the Macrobot system resulted in 158,469 data points. Each data point reflects one leaf measured. A total of 8,487 genotypes were phenotyped in 446 experiments. Data were analyzed in three hierarchical quality control steps. Two hundred ninety-nine of these data points were removed in the first step because they were considered outliers when comparing within technical replications. In the second step, 16 of the 446 experiments with a total of 5,571 data points were excluded because of the low infection rate of the susceptible control genotype (Kanzler) (Figure 2). In the third and final step, 697 data points were removed based on statistical outlier analysis among biological replications. Thus, 10,325 data points were excluded, representing only 6.5% of the total data. The overall data quality, as assessed by estimated heritability, was high. Data curation increased the heritability from 0.73 to 0.75 and reduced the variance component of the residuals (Table 1). We compared our data with results from the annual variety evaluation of the BSA. Macrobot data for PM resistance of 109 overlapping genotypes from the Elite panel showed a moderate correlation (r = 0.31). This correlation is within an expected range considering that our resistance assay is based on seedlings grown under greenhouse conditions and artificially inoculated with a specific PM isolate. In contrast, the BSA field trials are based on natural infections in field trials with locally adapted PM populations on adult plants.

The panel of 171 European elite wheat lines was on average more resistant than the 8,316 PGR. The mean percentage of infected leaf area of the elite lines and PGR were 25 and 48%, respectively. Nevertheless, 611 of the PGR were as resistant or even more resistant than the mean of the elite lines tested in our study (Figure 3). Furthermore, the Elite panel tested showed a binomial distribution of the BLUEs of infected leaf area. While our results illustrate the breeding progress in recent decades, the findings also highlight the need for the general implementation of PM resistance throughout the elite population and the potential of PGR as a source of this resistance.

The principal coordinate analysis (PCoA) performed using the GBS data revealed a narrow relationship among the 170 elite cultivars compared to the genetic diversity of the PGR panel (Figure 4) and the erosion of genetic diversity in the elite line pool. The distribution of the lines carrying beneficial alleles that may serve as donor lines is broad, reflecting the total genetic variation. Most of them show moderate or high genetic distance to the Elite panel. Although genetic diversity does not automatically imply a favorable phenotypic diversity, it demonstrates the potential of untapped genetic variation stored within the gene bank.

The GWAS performed based on the panel of wheat genetic resources resulted in 73 statistically significant marker-trait associations (MTA) located on 16 chromosomes (Figure 5; Tables 2–4; Supplementary Table 1). Two associated markers could not be assigned to a chromosomal location. The frequencies of the beneficial alleles (BAF) of the MTA were generally low in the population of PGR and absent or almost absent in the Elite panel. For the MTA, the mean frequency of the beneficial alleles was 0.004 in the Elite population. The ratios of explained variance (R²) of the 73 MTA ranged from 0.18 to 2% with mean R² of 1.1%. The 73 MTA were combined into 51 resistance loci based on flanking non-significant markers. These markers were in high LD within the defined regions (mean r² = 0.77). Individual loci spanned from 8 bp to almost 12 Mbp with a mean size of 0.8 Mbp, reflecting the wide variation in marker density of the subgenomes. Within all significant regions, 396 high confidence genes were predicted in the Wheat Reference Genome version 1.0 (Alaux et al., 2018). The number of genes per significant region varied from 0 to 50 (Supplementary Table 1), reflecting the differences in region sizes. For promising genes (Supplementary Table 1), we proposed KASP marker templates based on the reference genome for further use in research and breeding (Supplementary Table 6). To support the identification of candidate genes linked to markers significantly associated with mildew resistance, we performed domain enrichment analysis and observed highly significant enrichment of LRR domains in the identified significant regions as well as enrichment in other types of domains such as defensin, F-box-like domains among others (Supplementary Table 3). At this point, further detailed molecular genetic work is needed to reveal the genes and mechanisms behind the 51 discovered loci.

Identifying novel resistance loci in the PGR not currently used in European wheat breeding is one central research question. We compared the 51 MTA identified in this study with the extensive literature on reported PM resistance genes. We identified 11 MTA mapping to chromosome arms for which, to our knowledge, no PM resistance gene has been reported yet (Supplementary Table 4). These 11 MTA map to Chr3AL (MTA6-9), Chr7AS (MTA11), Chr3DL (MTA43-46) Chr4DL (MTA47), and Chr6DL (MTA49). The individual MTA differ significantly in size of the associated physical regions and the number of predicted genes located therein. Detailed information on these MTA and the putative candidate genes has been compiled in Supplementary Table 1. Several susceptible alleles were fixed in the commercial lines, leaving the breeding pool without an internal source of those favorable alleles. The PGR genotypes can provide the necessary resistance donors to compensate for the missing beneficial alleles in the elite line pool (Supplementary Table 2).

Host-pathogen interactions are a complex and dynamic process influenced by several biotic and abiotic factors. Naturally occurring infections can differ in population composition (Andrivon and De Vallavieille-Pope, 1993; Parks et al., 2009), and environmental factors (e.g., weather) play a critical role in their development (Cowger et al., 2012). Therefore, the use of experiment replications and the application of artificial infections can increase the precision in field phenotyping for resistance. Although this approach is feasible for smaller collections, its application to the entire PGR collection would require considerable efforts. Instead, the detached leaves assay combined with an automatic image analysis pipeline (Lück et al., 2020a) provided a high level of control over the environmental and biotic factors and ensured repeatable phenotypic data for the entire population at a reasonable cost in our study. The high observed heritability (0.75) support this expectation. Macrobot- and field data for PM resistance of 109 overlapping genotypes from the Elite panel showed a correlation of r = 0.31. Kang et al. (2020) observed a correlation of the PM resistance for 329 wheat varieties in the field and the greenhouse of 0.4–0.5. Tucker et al. (2006) reported an overall higher PM infection in the greenhouse than in the field using 293 recombinant inbred lines. In addition, when testing a BC1 double haploid population of 94 individuals, Mohler et al. (2013) observed correlations of 0.52 and 0.82 between adult plants in the greenhouse and in the field depending on the year. The latter study was based on a mixed pathogen population in the greenhouse and natural infection in the field. Moreover, high differences were also observed when comparing seedling and adult plant responses to PM (Wang et al., 2005; Jakobson et al., 2006; Mohler and Stadlmeier, 2019). Thus, we conclude that the method used in this study is less suitable for direct estimation of the genotype performance in the field, but rather provides the basis for the discovery of novel resistance mechanisms and crossing candidates for further trials. The Macrobot phenotyping provides a good proxy for fungal biomass (Lück et al., 2020b) and allows quantitative estimations of the fitness of a pathogen isolate on the host genotypes with only a small amount of seeds needed. Combined with a speed-breeding approach (Watson et al., 2018), this could significantly accelerate the breeding progress for resistance against PM. Another main difference to field conditions is the used pathotype. We used a single isolate, while the natural PM population structure is highly spatiotemporally diverse (Parks et al., 2009; Mascher et al., 2012) and responding to the deployment of resistance genes in the form of selection of favorable mutations and population shift (Andrivon and De Vallavieille-Pope, 1993). In this respect, the Macrobot system allows many genotypes to be tested against defined pathotypes, opening the opportunity to distinguish between monogenic R-gene resistance and durable race-nonspecific resistance mechanisms, thus assisting the integration of resistance alleles into elite varieties.

Powdery mildew resistance is a complex and long-studied trait with more than 90 genes at more than 60 loci reported so far (McIntosh et al., 2017; Tang et al., 2018; Ullah et al., 2018; Supplementary Table 4). However, only 12 of them have been cloned to date: Pm2 (Sánchez-Mart́ın et al., 2016), Pm3 (Yahiaoui et al., 2004), Pm4 (Sánchez-Mart́ın et al., 2021), Pm5 (Xie et al., 2020), Pm8 (Hurni et al., 2013), Pm17 (Singh et al., 2018), Pm21 (Cao et al., 2011; He et al., 2018; Xing et al., 2018), Pm24 (Lu et al., 2020), Pm38 (Krattinger et al., 2009), Pm41 (Li M. et al., 2020), Pm46 (Sánchez-Mart́ın et al., 2016), Pm60 (Zou et al., 2018),—reviewed also by Simeone et al. (2020). To demonstrate the potential value of the 51 resistance loci discovered in this study, we show a detailed examination of four of them (MTA5, MTA30, MTA37, and MTA50) which are interesting for different reasons. Briefly, MTA5 and MTA30 are co-localized with the well-known Pm4 and Pm26 locus, respectively. MTA37 is special because it is in contrast to the other MTA enriched in the tested Elite panel. MTA50 exhibit the highest estimated −log₁₀(p-value). The localization of the MTA was assumed to correspond to the Chinese Spring reference genome to enable comparison of the detected associations with known published genes. Of course, it is possible that the physical position of the MTA in the individual genotypes within the population differs from the reference genome, e.g., due to translocations, inversions, InDels, transposable elements, and other structural rearrangements of the genome (Appels et al., 2018; Bariah et al., 2020). Resistance genes in particular are subject to abundant structural variation, which decreases the mapping precision of associations and makes gene cloning thus more difficult (Dolatabadian et al., 2020; Nsabiyera et al., 2020). Most of the significant associations found were mapped to chromosome arms for which genes or QTLs associated with PM resistance have been previously reported. However, this does not necessarily mean that the associated marker is related to these genes. In many cases, the available sources do not include physical localization at a higher resolution than the chromosome arm or chromosome (Tables 2–4; Supplementary Table 1).

The many rare associations which are negatively fixed in the tested Elite panel show the huge potential in the tested PGR for modern wheat resistance breeding. We identified several accessions which carry multiple MTA from our experiment (see Supplementary Table 2). For all those MTA, the resistance conferring alleles are not present within the Elite panel and contribute significantly to the resistance against powdery mildew in our assay using the PGR. Those donors can be ordered from the Federal ex situ gene bank for Agricultural and Horticultural Crop Species in Gatersleben for prebreeding and scientific usage. To investigate the functionality in the field and use the described MTA in prebreeding we recommend the use of the Cadenza and Kronos TILLING Resources.

These and other studies have attempted to elucidate the genetic basis of resistance and isolate resistance genes for subsequent use in agriculture. However, with few exceptions, the identification and introgression of a single major QTL have not been confirmed to be a sustainable solution to ensure resistance as discussed previously in terms of the “boom-and-bust cycle” of resistance genes (McDonald and Linde, 2002). Furthermore, stacking many minor QTL is very laborious and time-consuming and even these quantitative resistances erode over time. Therefore, smarter strategies and interdisciplinary efforts are needed to protect our agroecosystems in the coming decades. One approach is to access novel resistance genes and mechanisms and optimize their management. McDonald (2014) proposed using dynamic diversity contrary to static diversity, based on deployment of single major-gene- or quantitative resistance, or constant R-gene pyramiding. However, this promising concept requires a significant increase in genotype diversity and breeding for populations composed of many genotypes. Diversifying the selection pressure exerted by the resistance genes used is one of the key strategies. As proposed by Hafeez et al. (2021), integrated approaches for creating a comprehensive wheat R-gene atlas should help cope with the challenge of the continuously changing pathogen populations. Another possibility would be establishing a classification in resistance genetics equivalent to the FRAC (Fungicide Resistance Action Committee) classification of fungicides. This classification groups all fungicides according to their mechanisms of action, allowing the farmer to diversify selection pressures. Growers could use such a classification to rotate and mix cultivars for resistance, as with fungicides. Regardless of the selected approach, the first and crucial step in the direction of resistance diversification is to reveal and exploit the genetic diversity stored in the gene banks worldwide. This information will allow growers and breeders to employ the new resistance genes even without knowing the mechanisms behind them.