Genomic signatures of selection for resistance to stripe rust in Austrian winter wheat

Key message We combined quantitative and population genetic methods to identify loci under selection for adult plant resistance to stripe rust in an Austrian winter wheat breeding population from 2008 to 2018. Abstract Resistance to stripe rust, a foliar disease caused by the fungus P. striiformis f. sp. tritici, in wheat (Triticum aestivum L.) is both qualitatively and quantitatively controlled. Resistance genes confer complete, race-specific resistance but are easily overcome by evolving pathogen populations, while quantitative resistance is controlled by many small- to medium-effect loci that provide incomplete yet more durable protection. Data on resistance loci can be applied in marker-assisted selection and genomic prediction frameworks. We employed genome-wide association to detect loci associated with stripe rust and selection testing to identify regions of the genome that underwent selection for stripe rust resistance in an Austrian winter wheat breeding program from 2008 to 2018. Genome-wide association mapping identified 150 resistance loci, 62 of which showed significant evidence of selection over time. The breeding population also demonstrated selection for resistance at the genome-wide level. Supplementary Information The online version contains supplementary material available at 10.1007/s00122-021-03882-3.


Introduction
Stripe rust is an economically important foliar disease of wheat (Triticum aestivum L.) caused by the fungus P. striiformis f. sp. tritici (Pst). Breeding resistant varieties are the most effective strategy for mitigating yield losses due to stripe rust (Chen 2020). Qualitative resistance to stripe rust in wheat is controlled by both qualitative resistance genes (R-genes) and quantitative trait loci (QTL) with small to moderate effects. More than 100 QTL have been associated with seedling resistance, adult plant resistance, and high temperature adult plant resistance in dozens of mapping populations and diversity panels (Rosewarne et al. 2008;Zegeye et al. 2014;Ye et al. 2019), and more than 80 Yr R-genes have been mapped or proposed to date (Waqar et al. 2018;Blake et al. 2019). Although Yr genes can provide complete or nearly complete protection against specific Pst races, they can easily break down with genetic changes in Pst populations (Poland et al. 2009;Chen 2020). For example, the emergence of the Warrior pathotype, which has overcome several widely deployed Yr genes, has caused devastating losses across Europe in the past decade (Buerstmayr et al. 2014;Hovmøller et al. 2016;Klymiuk et al. 2020;Tehseen et al. 2020). In contrast, small-and moderate-effect QTL provide partial, race non-specific resistance that tends to be more durable over time (Poland et al. 2009;Chen 2020). Information on Yr genes and QTL associated with stripe rust can be used for marker-assisted selection (Chen 2020) and to enhance genomic selection models (Juliana et al. 2017;Muleta et al. 2017).
Here, we analyzed historical stripe rust and genotyping data from an active Austrian winter wheat breeding program across 2008-2018. We employed genome-wide association (GWA) mapping to identify QTL associated with adult plant resistance to stripe rust within and across years and assessed their dynamics in allele frequencies and effects over the 11-year period to test for selection at the locus and genome-wide levels.

Phenotypic and genotypic data
Here, we analyzed a historical stripe rust dataset from the winter wheat breeding program of Saatzucht Donau GmbH & CoKG (Probstdorf, Austria). In total, 20,529 genotypes (12,844 recombinant inbred lines, 1638 doubled haploids, and 6047 advanced lines and registered varieties) were evaluated in 71 trials in 53 locations from 2008 to 2018 (Table 1). Because the plant material was part of an active breeding program, most genotypes were only evaluated in one plot in one trial (Table 1). To account for within-trial spatial variation, a check plot design was used, in which at least one genotype was replicated within each trial (Kempton 1984).
Each year, the Institute for Plant Protection in Field Crops and Grassland (Julius Kühn Institute, Kleinmachnow, Germany) provided urediniospores from a mixture of Pst pathotypes. The inoculum was then propagated on seedlings of susceptible genotypes in the greenhouse at the Saatzucht Donau research station in Reichersberg, Austria. One trial per year was grown in the Reichersberg disease nursery (location ID = LOC01), where plots were sprayinoculated with urediniospore suspension at the EC29 and EC30 growth stages (Leivermann and Brockerhoff 2015) (Online Resource 1). All other trials relied on natural infection. Each plot was scored for adult stripe rust resistance on a 1 (most resistant) to 9 (most susceptible) scale at 1-3 timepoints after symptoms became apparent on susceptible lines (Online Resource 1).
Genotypes (minimum F5 stage) with good agronomic performance (e.g., lodging resistance, yield, spike morphology), grain quality, and disease resistance (e.g., powdery mildew, Septoria nodorum blotch, stripe rust) were pre-selected for DNA sequencing. From the pre-selected material, a final subset of 5233 genotypes representing the diversity of the breeding program was chosen for sequencing and downstream genomic analysis (Table 1). Leaf samples from a minimum of ten plants per genotype were collected during early summer, and DNA was extracted as described by Saghai-Maroof et al. (1984). The DNA samples were genotyped with a custom 6 K Illumina marker array (Illumina, Inc., San Diego, CA, the USA) and with DArTseq (Diversity Arrays Technology Pty Ltd, Canberra, Australia) genotyping-by-sequencing (GBS) technology (Akbari et al. 2006;Elshire et al. 2011) and single nucleotide polymorphisms (SNPs) were then called using proprietary software. SNP genotypes were coded in terms of alternate alleles "a" and "A", where − 1 = aa (homozygous "a" allele), 0 = Aa (heterozygous), and 1 = AA (homozygous "A" allele). Missing SNP data was imputed with the "missForest" package (Stekhoven and Bühlmann 2012) in R (R Core Team 2020). To estimate imputation accuracy, 5% of the non-missing SNP data was masked (set to missing) and the dataset was imputed again, resulting in 94.4 ± 3.0% of correctly imputed masked SNPs. After filtering for minor allele frequency > 5% and call rate > 90%, a final set of 9744 SNPs was available for downstream genomic analyses (Online Resource 2). To generate a physical map, we used the nucleotide BLAST tool on the Wheat@URGI portal (Alaux et al. 2018) to compare the marker sequences against the IWGSC RefSeq v2.0 assembly (Appels et al. 2018). The physical position of each SNP was determined by the BLAST query with the greatest coverage value. 1 3

Phenotypic analysis
To adjust for spatial variation in each stripe rust score (stripe rust was scored at 1-3 different timepoints) within each trial, we fit a general linear model with genotype as a random effect and row and column effects modeled as two-dimensional P-splines and then estimated heritability using the "SpATS" package (Rodríguez-Álvarez et al. 2018) in R (R Core Team 2020). For the score with the greatest heritability in each trial, we fit a generalized linear model with genotype as a fixed effect and row and column effects modeled as two-dimensional P-splines with the "SpATS" package (Rodríguez-Álvarez et al. 2018) in R (R Core Team 2020) and then extracted the spatially adjusted stripe rust values (plot-level fitted values) for further analysis (Online Resource 1). We used the "lme4" package (Bates et al. 2015) in R (R Core Team 2020) to fit within-year (2013-2018) and across-year (2008-2018) mixed models with spatially adjusted stripe rust values as the response and genotype and trial as random effects. We extracted the variance components from each model and estimated broad-sense heritability (H 2 ) as ε is the error variance, and p h is defined as where n is the number of genotypes, and p i is the number of plots for the ith genotype (Holland et al. 2003). We also extracted the genotype best linear unbiased predictors (BLUPs) to estimate phenotypic correlations between years.

Genome-wise association
Because of the unbalanced nature of the dataset, we used methods that maximize statistical power for GWA in unbalanced studies (George and Cavanagh 2015;Xue et al. 2017;Chen et al. 2021). For within-year GWA, we used the onestage method, in which a mixed model is fit with plot-level phenotypes as the response and environmental (e.g., trial, year, location), genotypic (e.g., line, family), genetic background (e.g., relationship/kinship matrices, population structure components), and SNP information as fixed or random effects (Xue et al. 2017;Chen et al. 2021). Plant breeding experiments can include large numbers of individuals and/ or trials, making one-stage GWA computationally intensive when complex variance-covariance structures (e.g., relationship/kinship matrices) are included to control for background genetic effects (George and Cavanagh 2015;Xue et al. 2017). As such, the more common approach for GWA in plant systems has been the two-stage approach, in which (1) the plot phenotypes are regressed against environmental and genotypic terms and (2) the predicted genotypic means are then used as the phenotype in the GWA model including SNP and genetic background effects (George and Cavanagh 2015;Xue et al. 2017). Two-stage analysis can result in biased estimates in unbalanced studies, but methods have been developed to improve effect estimation when one-stage analysis is not computationally feasible (Möhring and Piepho 2009;Piepho et al. 2012;George and Cavanagh 2015;Xue et al. 2017). For across-year GWA, we employed a weighted two-stage analysis, which has been shown to closely approximate the results of one-stage analysis (Möhring and Piepho 2009;George and Cavanagh 2015;Xue et al. 2017).
For one-stage within-year (2013-2018) GWA, we fit mixed models with spatially adjusted stripe rust values as the response, SNP as a fixed effect, and genotype (only genotypes with SNP data) and trial as random effects using the "sommer" package (Covarrubias-Pazaran 2016) in R (R Core Team 2020). For within-year GWA from 2008 to 2012, the trial term was not included, as stripe rust was only evaluated in one inoculated trial in each of these years.
For two-stage across-year GWA, we first fit a mixed model with spatially adjusted stripe rust values as the response, genotype (all genotypes) as a fixed effect, and trial as a random effect using the "breedR" package (Muñoz and Sanchez 2020) in R (R Core Team 2020). We extracted the genotype best linear unbiased estimates (BLUEs) and standard errors (SE) of the genotype BLUEs from the model and then calculated the variances (σ 2 ) of the genotype BLUEs as where n is the number of observations per genotype (Online Resource 3). In the second stage, we used the "sommer" package (Covarrubias-Pazaran 2016) in R (R Core Team 2020) to fit a GWA mixed model with genotype BLUEs as the response, SNP as a fixed effect, and genotype (only genotypes with SNP data) as a random effect. For both within-year and across-year GWA, the variance of the genotype term was modeled as Kσ 2 a , where K is the realized additive relationship matrix (Endelman and Jannink 2012), and σ 2 a is the estimated additive genetic variance (Yu et al. 2006). For each model, we calculated K using SNP data from the genotyped lines included in the model with the "rrBLUP" (Endelman and Jannink 2012) package in R (R Core Team 2020). For across-year GWA, the residual variance was modeled as Iwσ 2 ε , where w is the vector of genotype BLUE variances (Möhring and Piepho 2009;George and Cavanagh 2015;Xue et al. 2017). The variance components were estimated once for each GWA model using the "population parameters previously determined" (P3D) method (Zhang et al. 2010).
Although K was included in all GWA models to account for population structure (Yu et al. 2006), there was little evidence of population structure in the breeding panel. We conducted a principal component analysis of the 5233 genotyped lines using SNP data with the "FactoMineR" (Lê et al. 2008) package in R (R Core Team 2020). The first and 1 3 second principal components accounted for 4.0% and 3.5% of the variance, respectively, and demonstrated some separation among lines with respect to the first year in which they appeared in the population (Online Resource 4).
SNP p values and effect estimates were extracted from each GWA model. For multiple test correction of the SNP p values, we conducted a false discovery rate (α = 0.05) analysis for each GWA model with the "qvalue" package (Storey 2015) in R (R Core Team 2020). The "sommer" package estimates SNP effect estimates (β) as = X � V − X X � V − y with X = ZM i , where Z is the incidence matrix of the genotype random effect, M i is the ith column of the SNP matrix, V − is the inverse of the phenotypic variance matrix, and y is the response (Covarrubias-Pazaran 2016). Because SNPs were coded as − 1 = aa, 0 = Aa, and 1 = AA, β is always relative to the number of "A" alleles.

Tests for selection
For each SNP, we calculated the frequency of allele "A" (p) in each year from 2008 to 2018 and extracted β from each within-year GWA model. To estimate the change in allele frequency of each SNP from 2008 to 2018, we fit a linear model for each SNP with p as the response and year as a continuous fixed effect and extracted the year coefficient from the model (Δp). Likewise, we estimated the change in allele effect of each SNP from 2008 to 2018 by fitting a linear model for each SNP with β as the response and year as a continuous fixed effect and extracted the year coefficient from the model (Δβ). In GWA, the power to detect a SNP-trait association and the absolute effect size of a SNP decrease with decreasing minor allele frequency (Bush and Moore 2012;Xiao et al. 2017). As such, effect sizes (1) become less negative (increase) as the major resistance allele increases in frequency and (2) become less positive (decrease) as the major susceptibility allele increases in frequency. To determine whether the frequency of the resistant allele or the susceptible allele of each SNP increased over time, we used the following criteria: (1) if Δp > 0 and Δβ < 0, there was selection for the "A" allele conferring susceptibility; (2) if Δp > 0 and Δβ > 0, there was selection for the "A" allele conferring resistance; (3) if Δp < 0 and Δβ < 0, there was selection for the alternate "a" allele conferring resistance; (4) if Δp < 0 and Δβ > 0, there was selection for the "a" allele conferring susceptibility.
We sought to test whether changes in allele frequencies were driven by selection rather than drift. For each SNP, we calculated the observed variance in allele frequency from 2008 to 2018 (V p ) and estimated the expected variance in allele frequency due to random genetic drift (V t ) as V t = p(1 − p) 1 − exp −t∕2N e , where p is the initial "A" allele frequency in 2008, t is the number of generations (t = 11 generations from 2008 to 2018), and N e is the effective population size (Ridley 2003;Juliana et al. 2019). We estimated N e (N e = 149) by regressing identity-bydescent (IBD) coefficients against time (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), with N e = 1∕2ΔIBD (Falconer and Mackay 1995). For each year (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), we calculated IBD between all pairs of lines using the SNPRelate package (Zheng et al. 2012) in R (R Core Team 2020). For each SNP, we then calculated the difference between the observed and expected variances, V p -V t . We compared V p -V t of each SNP to the genomewide null distribution of V p -V t . The null distribution was generated by subsampling V p -V t from 150 random SNPs in 1000 replications. The subsample size of 150 was selected because there were 150 significantly associated SNPs from GWA.
To test for genome-wide selection of stripe rust resistance or susceptibility, we estimated Ĝ , a composite statistic of the relationship between additive effect estimates and allele frequency changes over time of genome-wide markers, as described by Beissinger et al. (2018). We fit a random regression best linear unbiased prediction (rrBLUP) model with genotype BLUEs as the response (as described in the two-stage across-year GWA) and SNPs as fixed effects using the "rrBLUP" package (Endelman 2011) in R (R Core Team 2020). For each SNP, we extracted its estimated effect from the rrBLUP model and Δp (change in allele frequency from 2008 to 2018) from the selection analysis. We then estimated the value and significance of Ĝ with 1000 permutations using the "Ghat" package (Beissinger et al. 2018) in R (R Core Team 2020). As described by Beissinger et al. (2018), Ĝ = ∑ m j=1 Δ j j , where Δ j is the change in allele frequency from 2008 to 2018 for SNP j, α j is the allele effect of SNP j, and m is the total number of SNPs. To test whether the observed Ĝ was the result of selection rather than drift, Ĝ was compared to the null distribution of Ĝ perm (Beissinger et al. 2018). SNP allele effects were permuted 1000 times, and Ĝ perm was estimated for each permutation as Ĝ perm = ∑ m j=1 Δ j p j , where Δ j is the change in allele frequency from 2008 to 2018 for SNP j, p j is the allele effect of permuted SNP j, and m is the total number of SNPs (Beissinger et al. 2018). In this study, a negative Ĝ indicates selection for resistance to stripe rust and a positive Ĝ indicates selection for susceptibility.

Genotypic and trial effects on and heritability for stripe rust
From 2008 to 2012, stripe rust was evaluated on 962-1789 genotypes in one trial per year (Table 1). Stripe rust was evaluated on a larger panel of genotypes (1465-4134) in a greater number of trials (2-24) per year from 2013 to 2018 (Table 1). Broad-sense heritability (H 2 ) for resistance to stripe rust was generally high within years (H 2 = 0.50-0.90) and was moderate across years (H 2 = 0.54) ( Table 2). In most years, genotype explained a larger amount of the variance in stripe rust than trial and/or error (Table 2).
Between years, genotype BLUPs for stripe rust were positively correlated (Table 3). The number of genotypes in common was larger and phenotypic correlations tended to be stronger in adjacent years than in more distant years (Table 3). The highest correlations were observed between pairs of years from 2008 to 2012 (Table 3), where stripe rust was evaluated under artificial inoculation in the disease nursery in Reichersberg, Austria. From 2013 to 2018, trials were both artificially inoculated and naturally infected and were conducted in several locations. GWA across years and within 2009GWA across years and within -2011GWA across years and within , 2014GWA across years and within -2015GWA across years and within , and 2018 revealed 186 significant SNP-stripe rust associations (150 unique SNPs) after multiple test correction (Fig. 1, Online Resource 5-6). Of the significantly associated SNPs, 112 had a positive effect ("A" allele confers susceptibility) and 38 had a negative effect ("A" allele confers resistance) on stripe rust (Online Resource 6). The significant GWA SNPs explained a small proportion of the variance in stripe rust (R 2 = 0.08 ± 0.12) and had small to medium-effect sizes (|β|= 1.09 ± 1.23) (Online Resource 6).

Genome-wide association of stripe rust resistance
QTL colocalized between models at 12 locations (Fig. 1, Online Resource 6). The within-2010 and across-year GWA shared a SNP on chromosome 1A at 499.7 Mbp (Fig. 1, Online Resource 6). One SNP on chromosome 1D  identified (Online Resource 7). Few genotypes (N = 47) had SNP data in 2008 and the trial and residual terms explained a larger proportion of the variance in stripe rust in years with no significantly associated GWA SNPs, which may partially explain the lack of statistical power to detect SNP-trait associations (Online Resource 5).

Evidence of selection for stripe rust resistance
We assessed changes in allele frequencies and allele effects on stripe rust for each SNP from 2008 to 2018 and tested whether these changes were driven by selection or random genetic drift (Online Resource 8). By comparing the variance in observed allele frequencies to the expected variance due to drift (|V p -V t |) of each SNP against the null distribution of |V p -V t | (bootstrapped 1000 times, 95% quantile = 0.0008), we found significant evidence of selection of the resistant allele at 38/150 significant GWA SNPs ("A" allele at 23 SNPs; "a" allele at 15 SNPs) and for selection of the susceptible allele at 24/150 significant GWA SNPs ("A" allele at 8 SNPs; "a" allele at 16 SNPs) (Figs. 1, 2, Online Resource 8).
SNPs significantly associated in GWA and in selection tests demonstrated sharp changes in allele frequencies from 2012 to 2013, suggesting increased selection pressure during the generation between 2012 and 2013; the resistant or For SNPs with significant evidence of selection, the effect (A) and frequency (B) of the allele under selection (regardless of "A" or "a" allele state) are plotted against time, with SNPs with selection for the resistant allele in blue and for the susceptible allele in red. For SNPs not under selection, the effect (C) and frequency (D) of the major allele (allele at higher frequency, regardless of "A" or "a" allele state) are plotted against time, with blue and red denoting resistance and susceptibility conferred by the major allele, respectively 1 3 susceptible allele of these SNPs was nearly fixed in the population by 2018 (Fig. 2, Online Resource 8). In contrast, the allele frequencies and effects of the significant GWA SNPs that were not under selection were relatively unchanged from 2008 to 2018. Of the 88 significant GWA SNPs not under selection, the major allele (allele at higher frequency, regardless of "A" or "a" allele state) conferred resistance at 71 SNPs and susceptibility at 17 SNPs (Fig. 2, Online Resource 8). Allele effect estimates may have been inflated in 2008-2012, as stripe rust was only evaluated in one trial per year and fewer genotypes had SNP data in these years than in 2013-2018 (Fig. 2).
To test whether genome-wide resistance or susceptibility to stripe rust was under selection in the breeding program between 2008 and 2018, we used the Ĝ method (Beissinger et al. 2018). There was significant evidence of genomewide selection for stripe rust resistance over the 11 years in the population, as demonstrated by a negative Ĝ value ( Ĝ = − 0.26) and a highly significant (p = 2 × 10 -16 ) difference between the observed Ĝ and the null distribution of 1000 permuted Ĝ perm values (Fig. 3A).
SNPs with larger effect sizes on stripe rust in across-year GWA had greater changes in allele frequencies from 2008 to 2018 (Fig. 3B). Furthermore, 2483 SNPs were not significantly associated with stripe rust in GWA within or across years, yet they had significant evidence of selection for the resistant (1233 SNPs) or the susceptible (1250 SNPs) allele over time (Online Resource 8). For SNPs under selection, absolute allele effect sizes from GWA within and across years (|β|) and V p -V t were greater at significant GWA SNPs (|β|= 0.42 ± 0.36; V p -V t = 0.07 ± 0.06) than at nonsignificant GWA SNPs (|β|= 0.12 ± 0.15; V p -V t = 0.005 ± 0.006). These results suggest that both moderate and small effect QTL were under selection, although to a lesser extent for QTL with small effects that were not detectable GWA.

Discussion
Selection pressure for stripe rust resistance can be influenced by both breeder's decisions and changes in pathotype composition of Pst populations. Here, we combined quantitative genetic and population genetic methods to identify genomic regions that were under selection for resistance or susceptibility to stripe rust in an Austrian winter wheat breeding program from 2008 to 2018. GWA revealed 150 SNPs significantly associated with stripe rust within 2009-2011, 2014-2015, and 2018 and across 2008-2018, many of which overlapped with regions previously associated with stripe rust resistance in other populations (Rosewarne et al. 2013) and with putative Yr R-genes (Waqar et al. 2018;Blake et al. 2019). Because the ability to detect SNP-trait associations is largely dependent on minor allele frequency (Bush and Moore 2012;Xiao et al. 2017) and selection within a breeding program can generate rapid changes in allele frequencies (Ridley 2003), the majority of these SNPs were detected by GWA in only 1 year or in adjacent years. Investigating the dynamics in allele frequencies and effects over time can identify regions of the genome which have undergone selection for specific traits (Juliana et al. 2019). By combining GWA and selection testing, we found that both small-and moderate-effect loci had evidence of selection in the population. We also employed the Ĝ method to assess selection at the genome-wide level (Beissinger et al. 2018) and found that the breeding population demonstrated genome-wide selection for resistance from 2008 to 2018.
The highly significant QTL on the short arm of chromosome 2A were under selection for resistance and although it was physically near the Yr17 gene (Rosewarne et al. 2013), it is unlikely that Yr17 underlies this QTL because virulent Pst races overcame Yr17 across European wheat cultivars by 2000 (Bayles et al. 2000). Two significant GWA SNPs on the short arm of chromosome 2B also demonstrated selection for the resistant allele and may be linked to Yr27, an R-gene which has recently broken down against new Warrior-type races of Pst in the Middle East (Tehseen et al. 2020). A SNP in the pericentromeric region of chromosome 1A was under selection for the resistant allele and was near QTL for adult plant resistance to stripe rust from four mapping populations (Rosewarne et al. 2008;Dedryver et al. 2009;Bariana et al. 2010;Ren et al. 2012) and in a panel of elite spring wheat lines from CIMMYT (Crossa et al. 2007), but no Yr genes have been mapped to this region (Waqar et al. 2018;Blake et al. 2019). Four SNPs on the long arm of chromosome 5B were under selection for susceptibility and colocalized with QTL associated with non-race-specific adult plant resistance to stripe rust in a Sichuan wheat diversity panel (Ye et al. 2019) and with QTL for race-specific seedling resistance to stripe rust found in two biparental mapping populations (Feng et al. 2011;Zegeye et al. 2014) and for adult plant resistance in an Austrian biparental mapping population (Buerstmayr et al. 2014). However, no Yr genes have been mapped to the long arm of chromosome 5B to date (Waqar et al. 2018;Blake et al. 2019). Two SNPs on the long arm of chromosome 7A were also under selection for susceptibility, but we found no evidence of previously reported stripe rust QTL or mapped Yr genes in this region (Waqar et al. 2018;Blake et al. 2019).
By combining SNP-specific and genome-wide approaches, we demonstrated that the breeding population harbors both moderate-effect QTL and quantitative forms of incomplete, race non-specific adult plant resistance and that both were under selection across the 11-year period. The resistance QTL identified in this study will be further evaluated for their use in marker-assisted selection and as covariates in genomic prediction models for stripe rust resistance in the breeding program. The breeding population demonstrated highly heritable, quantitative resistance to stripe rust and low population structure, indicating that genomic prediction of stripe rust resistance can be successfully applied in this population (Crossa et al. 2017;Juliana et al. 2017;Muleta et al. 2017;Tehseen et al. 2021).
Author Contribution statement LM, SM, HB contributed to conceptualization; LM, SM, CA were involved in data curation; LM, SM, HGD, SZ contributed to formal analysis; HB was involved in funding acquisition; FL, AN contributed to investigation; LM was involved in methodology, validation, visualization, writing-original draft; LM, SM, FL, AN, HB contributed to project administration; FL, AN, CA, HB were involved in resources; LM, SM, HGD, SZ, CA contributed to software; SM, FL, HB were involved in supervision; LM, SM, HB contributed to writing-review and editing. Availability of data and materials All phenotypic and genotypic data and results from the analyses presented here are included in the manuscript materials.

Code availability
The scripts used to conduct the analyses presented here are available upon request.
Conflict of interest FL, AN, and CA were employed by the company Saatzucht Donau GmbH & CoKG. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.