Separation of the effects of two reduced height (Rht) genes and genomic background to select for less Fusarium head blight of short-strawed winter wheat (Triticum aestivum L.) varieties

Key message FHB resistance shared pleiotropic loci with plant height and anther retention. Genomic prediction allows to select for genomic background reducing FHB susceptibility in the presence of the dwarfing allele Rht-D1b. Abstract With the high interest for semi-dwarf cultivars in wheat, finding locally adapted resistance sources against Fusarium head blight (FHB) and FHB-neutral reduced height (Rht) genes is of utmost relevance. In this study, 401 genotypes of European origin without/with dwarfing alleles of Rht-D1 and/or Rht24 were analysed across five environments on FHB severity and the morphological traits such as plant height (PH), anther retention (AR), number of spikelets per ear, ear length and ear density. Data were analysed by combined correlation and path analyses, association mapping and coupling single- and multi-trait genome-wide association studies (ST-GWAS and MT-GWAS, respectively) and genomic prediction (GP). All FHB data were corrected for flowering date or heading stage. High genotypic correlation (rg = 0.74) and direct path effect (0.57) were detected between FHB severity and anther retention (AR). Moderate correlation (rg = − 0.55) was found between FHB severity and plant height (PH) with a high indirect path via AR (− 0.31). Indirect selection for FHB resistance should concentrate on AR and PH. ST-GWAS identified 25 quantitative trait loci (QTL) for FHB severity, PH and AR, while MT-GWAS detected six QTL across chromosomes 2A, 4D, 5A, 6B and 7B conveying pleiotropic effects on the traits. Rht-D1b was associated with high AR and FHB susceptibility. Our study identified a promising positively acting pleiotropic QTL on chromosome 7B which can be utilized to improve FHB resistance while reducing PH and AR. Rht-D1b genotypes having a high resistance genomic background exhibited lower FHB severity and AR. The use of GP for estimating the genomic background was more effective than selection of GWAS-detected markers. We demonstrated that GP has a great potential and should be exploited by selecting for semi-dwarf winter wheat genotypes with higher FHB resistance due to their genomic background. Supplementary Information The online version contains supplementary material available at 10.1007/s00122-022-04219-4.


Introduction
Wheat (Triticum aestivum L.) is one of the most cultivated cereal crops worldwide and serves as a staple food for millions of people. However, the production of wheat is hampered by several diseases, including Fusarium head blight (FHB) caused by numerous Fusarium species, with Fusarium graminearum Schwabe and F. culmorum (W.G. Smith) Sacc. being the most predominant in Central Europe. FHB causes severe yield losses and contaminates the grains with several mycotoxins including deoxynivalenol (DON) which is one of the most frequently detected in wheat (Righetti et al. 2021;Topi et al. 2021). DON makes the grain unsuitable for flour and malt and is also toxic for non-ruminant animals (Windels 2000). Damages due to FHB in wheat are likely to increase with the rising temperatures and higher atmospheric carbon dioxide (CO 2 ) caused by climate change (Timmusk et al. 2020;Miedaner and Juroszek 2021;Hay et al. 2022). Fungicides and DON-reducing technologies used as traditional measures to control FHB disease increase production costs, with no significant positive return on grain yield (Wilson et al. 2018). Most effective Communicated by Hermann Buerstmayr. 1 3 and environmentally friendly strategies integrate appropriate agronomic practices and enhanced host resistance.
Breeding for resistance against FHB in European winter wheat faces several challenges related to the complex genetic architecture of the trait (Jiang et al. 2015;Ruan et al. 2020) and requires intensive breeding effort and accurate breeding strategies . FHB resistance is controlled by multiple medium-and small-effect quantitative trait loci (QTL) which are vulnerable to the changing environmental conditions due to QTL-environment interactions (Miedaner and Juroszek 2021). To date, more than 500 QTL have been reported for FHB from different breeding populations and were further clustered into 65 QTL with unique physical positions through a recent meta-analysis by Venske et al. (2019). Additionally, morphological traits such as plant height and anther retention or extrusion were shown to passively contribute to FHB resistance (Mesterházy 1995;Buerstmayr et al. 2020;Xu et al. 2020), and relationships between FHB resistance and morphological traits have been extensively reported using correlation coefficients (Buerstmayr and Buerstmayr 2015;Steiner et al. 2019b;Ruan et al. 2020;Xu et al. 2020). Although correlation analysis helps measure the degree of association between two traits, it does not explain causes and effects of the relationship (Dewey and Lu 1959;Ozukum et al. 2019). In addition, the existence of strong correlation between two traits might be due to the presence of one or many other traits which strengthen the complexity of the interactions. Understanding the complex interactions between FHB resistance and morphological traits requires the use of appropriate statistical approaches like path analysis, also referred to as structural equation modelling (SEM) or covariance structural equation modelling (CSEM).
Path analysis is a causal multivariate modelling approach which complements simple correlation analysis by unravelling the nature and magnitude of the observed relationships among traits (Wright 1934). It exploits observed correlations to estimate standardized direct and indirect effects contrary to the standard multivariate modelling which ignores causal relationships among variables and combines all effects together (Valente et al. 2013). Path effects estimated using correlations are unitless and interpreted as standardized regression slopes, allowing for comparison of the relative importance of different variables (Stage et al. 2010). Therefore, a direct path effect of a particular morphological trait on FHB resistance would indicate how much an increase of one unit in the standard deviation of that trait directly changes the standard deviation of FHB resistance. On the contrary, an indirect path effect would indicate how much a particular morphological trait changes the standard deviation of FHB resistance depending on the presence of other morphological traits. The selection for FHB resistance using the indirect path of a trait requires the simultaneous integration of one or several other traits. Based on comparison between the two types of effect, the breeder can decide which paths and morphological traits are more important to be considered for improving FHB resistance using correlated traits for higher genetic progress. Hence, a combined correlation and path analysis has several advantages including setting reliable criteria for multiple-trait selection, minimizing risks of components compensation and guiding in planning of experiments (Usman et al. 2017;Khan et al. 2022). Path analysis was used by researchers to depict complex associations among yield and related traits in several crop species including wheat (Baye et al. 2020;Hinson et al. 2022) and maize (Toebe and Cargnelutti 2013).
Plant height in wheat is controlled by at least 150 QTL scattered across the whole wheat genome (Mao et al. 2010), and 25 Rht genes with more than 30 dwarfing alleles are reported (Sanchez-Garcia and Bentley 2019). They became very popular after their introduction during 1960s and 1970s to initiate the "Green Revolution" in developing countries (Hedden 2003) but also to breed short-strawed wheat in industrial countries. Major Rht genes are Rht-B1 (syn. Rht1), Rht-D1 (syn. Rht2) and Rht24, located on chromosomes 4B, 4D and 6A, respectively (Würschum et al. 2017;Tian et al. 2022). The wild type is named "a" and the height-reducing mutant "b". Rht-D1b and Rht24b were more frequently found in Central European wheat germplasm than Rht-B1b. Rht genes have been widely used in commercial wheat breeding to develop semi-dwarf varieties which are preferred to tall genotypes because of their lodging resistance, higher nitrogen fertilizer use efficiency, increased tillers, higher harvest index and grain yield (Zhao et al. 2022).
Unfortunately, Rht-B1b and Rht-D1b were associated with FHB susceptibility and lower anther extrusion (Lu et al. 2013;Hales et al. 2020). Miedaner et al. (2019) and Zhang et al. (2021) suggested the use of major FHB-resistant QTL through marker-assisted introgression to counterbalance the negative FHB effect of semi-dwarf genotypes. Meanwhile, it has also been demonstrated that major non-adapted FHB resistance QTL (Fhb1, Fhb2, Fhb5) provide only a partial resistance (Su et al. 2019), and their introgression was not equally effective across different backgrounds (Brar et al. 2019) and requires considerable breeding efforts. Breeders must also consider a larger number of medium-and smalleffect QTL in locally adapted cultivars for all traits analysed here, so-called genomic background (Brar et al. 2019) in order to aim for local, short-strawed wheat cultivars with higher FHB resistance. The genomic background refers to all locally adapted genes which are likely to influence the phenotype controlled by a particular gene of interest (Yoshiki and Moriwaki 2006). In wheat, all genes of a genotype which are capable of influencing the reduction effect of major Rht genes on plant height can be considered as genomic background for the plant height in that genotype and similarly for FHB resistance and AR. Numerous medium-and smalleffect loci have been reported for FHB resistance and anther retention, but their potential to counterbalance the negative effect major Rht genes such as Rht-D1 and possible exploitation in practical breeding have been poorly investigated. Herter et al. (2018) evaluated segregating materials derived from a biparental cross between two German winter wheat cultivars, "Solitär x Bussard", and demonstrated that Rht24b was not associated with FHB resistance. This was recently confirmed in a large European winter wheat germplasm by Miedaner et al. (2022) who also found that Rht24b did not affect FHB resistance. Therefore, strategies to develop FHB-resistant varieties with semi-dwarfing alleles should reconsider the choice of such Rht genes. This gives the opportunity to develop marker arrays and strategies to facilitate the exploitation of FHB-neutral Rht genes in breeding programmes. However, given the complex interactions among morphological traits and their passive contribution to FHB resistance ), a clear understanding of the effect of Rht24b on morphological traits, like anther retention or extrusion, is required. Furthermore, the choice of Rht genes could also be driven by breeding objectives and the extent to which FHB-neutral Rht genes reduces the plant height compared with other genes.
In this study, we aimed to (i) describe the nature and magnitude of the complex interactions between FHB resistance and morphological traits, including plant height, anther retention, number of spikelets per ear, ear length and ear density using combined correlation and path analyses; (ii) dissect the genetic architecture of FHB severity, plant height and anther retention, thus uncovering the genetic basis of the complex interactions among these traits through a joint implementation of single-and multi-trait genome-wide association study (GWAS), and genomic prediction (GP); and (iii) separate for each trait, the effects of Rht genes and genomic background (GB), evaluating the potential of GB in selecting FHB-resistant genotypes with short Rht alleles based on the single-trait GWAS and genetic estimated breeding values (GEBV) from the GP. GB refers to all QTL affecting the traits, except those associated with plant height on chromosomes 4D and 6A corresponding to Rht-D1 and Rht24 genes, respectively.

Plant materials and field experiments
The materials consisted of 401 winter wheat cultivars from European origin (Table S1). These cultivars were evaluated in Germany at three locations in 2020 and two in 2021, resulting in five environments (location × year combinations) in total. In 2020, the cultivars were evaluated in Hohenheim (HOH) near Stuttgart (9.12° E, 48.42° N; 400 m above sea level [a.s.l]), Oberer Lindenhof (OLI) near Reutlingen (9.18° E, 48.28° N; 700 m a.s.l), and Wohlde (WOH) near Bergen (9.98° E, 52.80° N; 80 m a.s.l). In 2021, field experiments were conducted in HOH and OLI only. At each location, the cultivars were randomized using an incomplete lattice design with two replicates. Experimental units were planted in double rows of 0.9 m length in WOH and single rows of 1.2 m length in HOH and OLI. Experiments were sown with a density of 40 kernels/row in WOH and 60 kernels/row in HOH and OLI.

Artificial inoculations
Inoculum of the highly aggressive single-spore isolate FC46 (IPO 39-01) of F. culmorum (Snijders and Perkowski 1990) was produced on autoclaved wheat kernels as described in detail by Miedaner et al. (1995Miedaner et al. ( , 1996. Prior to inoculation, the Fusarium suspension was diluted to a concentration of 2 × 10 5 spores mL −1 . Approximately 100 mL m −2 of the diluted inoculum was applied using an adapted agricultural sprayer (Hege 75, Waldenbuch, Germany). Inoculations were repeated four to five times at intervals of two to three days to inoculate each genotype at least once during midanthesis. The first inoculations were done when early cultivars started flowering.

Phenotypic data collection
Eight traits were recorded: Fusarium head blight severity (FHB severity, %), plant height (PH, cm), anther retention (AR, %), ear length (EL, cm), number of spikelets per ear (NS, spikelets ear −1 ), ear density (ED, spikelets cm −1 ), days to flowering (DTF, days) and heading stage (HS). Days to flowering, the number of days when 75% of heads showed extruded anthers after May 1, was assessed in all environments except in WOH. Heading stage (1-9) was recorded on one a single day (different dates in different environments) and is equal to BBCH stage 51-59 (Meier 2001). Plant height was recorded for each plot after flowering. FHB severity was rated as the percentage of infected spikelets (0-100%) of each head, averaged across all plants of a plot. The first rating started at the onset of symptom development about 15-20 days after inoculation and was repeated in intervals of three to five days until the first signs of ripening. We rated each time the total number of infected spikelets per plot. This rating approach at different stages in the development of the plants under field conditions helped to evaluate the combination of both type I (incidence) and type II (symptom development) FHB resistance in one number (Nannuru et al. 2022). Anther retention was assessed as percentage of retained anthers per spike according to Atashi-Rang and Lucken (1978). Five spikes per plot were harvested five to 1 3 seven days after flowering and frozen in paper bags to do the assessment in off-season. Retained anthers of two lateral florets of eight central spikelets were counted on each of the five spikes. In contrast to Atashi-Rang and Lucken (1978), anthers located between the lemma and palea, and partially extruded, were counted as retained. From the same spikes, the number of spikelets per ear (NS i ) was counted and ear (spike) length (EL i , cm) was measured. Ear density ED i was calculated as the number of spikelets per centimetre: AR, EL and ED were estimated plot-wise by averaging values over the five spikes.

Genotyping and molecular analysis
DNA isolation and genotyping were done by SGS Institut Fresenius, TraitGenetics Section (https:// trait genet ics. com/ index. php/ conta ct, Gatersleben, Germany) from seed samples of the 401 genotypes. Genotyping was done using an Illumina 25 K Infinium single-nucleotide polymorphism (SNP) array which produced a total of 24,145 SNP markers. About 58 and 14% of the markers overlapped, respectively, with 90 K iSelect and 820 K Axiom ® arrays, which are publicly available at CerealsDB (http:// www. cerea lsdb. uk. net/ cerea lgeno mics/ Cerea lsDB) (Wilkinson et al. 2020). SNPs were filtered by removing markers with minor allele frequency (MAF) < 0.05 and call rate < 80%. This narrowed down the marker data to 19,969 polymorphic and high-quality SNPs with a total of 0.7% of missing data. SNPs were coded as − 1, 0 and 1 corresponding to AA, Aa and aa. A represented the major allele, while a denoted the minor allele. Heterozygous loci were replaced with missing values, and the data were imputed using the Wright's equilibrium approach (Wright 1922). SNPs filtering, coding and imputation were conducted using the raw.data function in the snpReady R package (Granato et al. 2018). Start and end physical positions (Table S6) of the markers sequences were obtained by blasting against the International Wheat Genome Sequencing Consortium (IWGSC) reference sequence (IWGSC RefSeq) v.2.1 which is accessible at the Research Unit in Genomics-Info (URGI) of the French National Institute for Agriculture, Food and Environment (INRAE) (https:// wheat-urgi. versa illes. inra. fr/ Seq-Repos itory/ Assem blies) .

Phenotypic data analysis
For each plot, arithmetic mean of repeated ratings of FHB severity was calculated and included in the statistical analysis as described by Mesterházy (1995). Estimated individual value was validated when the corresponding coefficient of variation of error was below 5%. Outliers were removed from the data using the Bonferroni-Holm method (Bernal-Vasquez et al. 2016). Variance components, genetic coefficient of variation, broad-sense heritability and best linear unbiased estimations (BLUEs) were calculated for each of the eight traits, based on a mixed linear model: where y ijkl is the phenotype of the ith genotype in the jth environment within the kth block of the lth replicate; g i is the effect of the ith genotype; e j is the effect of the jth environment; ge ij is the effect of the interaction between the genotype and the jth environment; b k is the effect of the kth block; r l is the effect of the lth replicate; and ε ijkl is the residual error on the phenotype of the ith genotype in the jth environment within the kth block of the lth replicate. To adjust for the impact of repeated artificial inoculations on early flowering genotypes (higher inoculum dose after flowering), days to flowering was added as a cofactor (fixed effect) to the mixed linear model for calculating FHB severity. Because days to flowering was not assessed in WOH, the model was further extended by adding heading stage as a second cofactor in order to correct FHB severity. This correction reduced the correlation coefficient between FHB severity and days to flowering from − 0.32 to − 0.16. Similarly, the correlation between FHB severity and heading stage was reduced from 0.29 to 0.11. The genotype was treated as fixed effect, while block, replicate, environment and genotype-environment interaction were included as random effects. Variance components were firstly estimated for all genotypes as one group. Secondly, genotypes were grouped based on Rht genes and variances components were estimated to describe the extend of genetic variation within each group for PH, FHB severity and AR. Prior to genotypes grouping, a singletrait genome-wide association study was conducted on each trait to identify significant markers associated with PH on chromosomes 4D and 6A, which correspond to Rht-D1 and Rht24 genes. Based on the presence/absence of Rht-D1b and Rht24b, genotypes were categorized into four Rht groups, namely NoRht, Rht24b, Rht-D1b and Rht24b + Rht-D1b. NoRht was composed of genotypes with tall alleles of both genes. Rht24b included genotypes carrying Rht24b only, while Rht-D1b was constituted of genotypes with Rht-D1b only. Similarly, Rht24b + Rht-D1b included genotypes possessing semi-dwarfing alleles of both genes. Dummy variables (0, 1) as described by Piepho et al. (2006) were created to separate genotypes, and genotype-dummy interaction was added to the model to estimate variance components within each Rht group. For dummy = 1, the interaction between genotype and dummy estimates variance components within the respective group. Significance of variance components (2) y ijkl = g i + e j + b k + r l + ge ij + ijkl 1 3 was calculated using a likelihood ratio test (LRT) based on factor-wise comparisons of a model with all factors and a model without the respective factor (Stram and Lee 1994;Morrell 1998). Environment-specific (heterogeneous) variances were fitted for replicate, block and residual effects. Harmonic means of environment-specific residual variances were calculated and reported. The mixed linear models were fitted using the ASReml-R package v.4.1.0 (Gilmour et al. 2015). Only adjusted means (BLUEs) of the corrected FHB severity and other traits across the five environments (location × year combinations) were considered for further analyses (Table S1).
Broad-sense heritability (H 2 ) was estimated for each trait following the procedure described by Piepho and Möhring (2007): where 2 G is the genotypic variance and v BLUE Δ… is the mean variance of the difference of two adjusted means (BLUEs). The genetic coefficient of variation (CV G ) was calculated for each trait as follows: where 2 G is the genotypic variance and x is the mean of BLUEs of the trait. Coefficient of variation of error (CV ε ) was also calculated for each trait by replacing the genotypic variance by the residual variance ( 2 ) in Eq. (4). All analyses were performed in R software, v.4.1.3 (R Core Team 2021).

Correlation and path analysis
Genotypic correlations were calculated between pairs of traits by fitting a bivariate mixed linear model as follows: where y ijkl is the phenotypic value of the first trait; y' ijkl is the phenotypic value of the second trait; g i is the effect of the ith genotype; e j is the effect of the jth environment; ge ij is the effect of the interaction between the ith genotype and the jth environment; b k is the effect of the kth block; r l is the effect of the lth replicate; and ε ijkl is the residual error. Bivariate models were fitted with unstructured variance-covariance using the corh-option for the genotype, environment and genotype by environment interaction and diagonal variance-covariance for replicate and blocks. In addition to genotypic correlations, Pearson's productmoment correlations analysis was performed based BLUEs and the results plotted using the GGally package v.2.1.2. Correlation coefficients were classified and interpreted as described by Zou et al. (2003).
The genotypic correlations were used to perform path analysis between FHB severity, as response variable, and morphological traits such as PH, AR, EL and NS included as explanatory variables in agricolae R package v.1.3-5 (De Mendiburu 2016). The analysis was based on standardized partial regressions which used observed genotypic correlations to estimate direct and indirect path effects as illustrated diagrammatically in detail in Fig. 1 (Wright 1918(Wright , 1934Dewey and Lu 1959). For simplicity, PH, AR, EL and NS were coded as 1, 2, 3, and 4, respectively. The path model was fitted as follows: where ̂ is the standardized regression slope, representing the direct effect estimators; Y is the vector of genotypic correlations between FHB severity and morphological traits; X is the genotypic correlation matrix among morphological traits; and ε is a vector of residual errors. To estimate the direct effect of each exploratory variable on FHB severity, we minimized the residual least squares and took its derivative with respect to ̂ , creating a system of normal equations as described by Toebe and Cargnelutti (2013): where X' and X −1 are the transpose and inverse of the correlation matrix X, respectively; p 15 , p 25 , p 35 and p 45 are the direct path effects of PH, AR, EL and NS on FHB severity, respectively. Afterwards, the indirect path effects were estimated using structured equations (Dewey and Lu 1959) linking the correlation coefficients and the direct path effects as follows: where r i5 is the genotypic correlation coefficient between the ith trait and FHB severity; p i5 is the direct path effect of the ith trait on FHB severity; and ∑r ik p k5 is the summation of indirect path effects of the ith trait on FHB severity via all other k traits. We obtained four structured and simultaneous equations which defined a matrix of direct (diagonal elements) and indirect path (off-diagonal elements) effects of morphological traits on FHB severity (Fig. 1b). Coefficient of determination (R 2 ) of the path model was estimated as: where r 15 , r 25 , r 35 and r 45 are the genotypic correlations between FHB severity and PH, AR, EL and NS, respectively. The direct effect (p res ) of the residual factors was calculated as the part of the variation of FHB severity which was not explained by PH, AR, EL and NS (Toebe and Cargnelutti 2013). That is, The residual factors referred to all other variables which contribute to the variation of FHB severity and were not included in the path model. Similar to the standard multivariate modelling, the existence of high multicollinearity among exploratory variables can lead to significant biases in the estimates of path effects and result in a wrong interpretation of existing causal relationships among the variables (Stage et al. 2010;Toebe and Cargnelutti 2013). Therefore, we investigated the multicollinearity among exploratory variables by calculating the variance inflation factor (VIF). In doing this, a linear regression was fitted for each variable on other variables, and the respective multiple R 2 ( R 2 mult ) was used to calculate VIF as: Multicollinearity was considered as high, if VIF > 5 (Olivoto et al. 2017).

Population structure and genome-wide association studies
Polymorphic and high-quality SNP markers obtained after the filtering were used to investigate the population structure of the 401 genotypes, using the snpdgsPCA() function available in SNPRelate R package v.1.26 (Zheng 2015). Scatter plots illustrating the structure of the population were built for the first two principal coordinates using the ggplot2 package v.3.3.5 (Wickham et al. 2021).
To dissect the genetic architecture of FHB severity, PH and AR, single genome-wide association study (ST-GWAS) analysis was performed based on all 19,969 SNP markers with MAF ≥ 0.05 to identify significant marker-traits associations (MTAs), using the Bayesianinformation and Linkage-disequilibrium Iteratively Nested Keyway (BLINK) method available in the Genomic Association and Prediction Integrated Tool (GAPIT) R package v.3.1.0 (Wang and Zhang 2021). BLINK is a multi-locus GWAS method which was demonstrated to outperform in both speed and statistical power other methods such as the fixed and random model circulating probability unification (FarmCPU) and mixed linear model (MLM) (Huang et al. 2019;Wang and Zhang 2021). In contrary to FarmCPU, BLINK works directly on the markers instead of bins, thereby removing the assumption that causal loci are evenly distributed across the genome. BLINK includes two fixed effects models (FEM 1 and FEM 2 ) and one filtering process, which are performed iteratively (Huang et al. 2019). The filtering process consists in sorting all markers and selecting the first t most significant markers referred to as pseudo-quantitative trait nucleotides (QTNs) (based on indicated threshold) that are not in linkage disequilibrium (Pearson's correlation < 0.7) as covariates. Based on this, the two FEMs can be written as follows: where y i is the BLUE across environments of the ith genotype; S i1 , S i2 , …, S ik are the genotype scores of k pseudo-QTNs, initiated as an empty set; b 1 , b 2 , …, b k are the respective effects of the pseudo-QTNs; S ij is the genotype score of the ith genotype and jth marker; d j is the corresponding effect of the jth marker; and e i is the residual with e i ~ N(0, 2 ). In FEM 1 , S ij d j represents a testing markers term where all remaining markers are tested one at a time. In FEM 2 , all t pseudo-QTNs selected after the filtering are re-examined together using the Bayesian information content (BIC) as an optimization criterion. Afterwards, k out of the t pseudo-QTNs with the best BIC values are selected and included in FEM 1 as covariates to test the remaining markers.
As BLINK method does not automatically make a clear separation between major-and small-effect QTL (pseudo-QTNs selected are retested together and some are also removed in FEM 2 ), we implemented two ST-GWAS analyses (ST-GWAS 1 and ST-GWAS 2 ) on each of the trait. In ST-GWAS 1 , all markers were included in the analysis, while in ST-GWAS 2 significant MTAs identified on chromosomes 4D (Rht-D1) and 6A (Rht24) from ST-GWAS 1 were used as covariate variable to increase the detection of small-effect markers. In wheat, the use of Rht markers as covariate variables in BLINK was reported to increase the number of significant markers (Merrick et al. 2022). Significant MTAs were identified based an exploratory threshold of − log10 (p value) = 6 and a Bonferroni-corrected threshold of − log10 (p value) = 6.30. The corrected Bonferroni threshold (p value) was determined as follows: where α = 0.01 is the type 1 error and Me is the number of markers included in the analysis. The total proportion of genetic variation (p G ) explained by significant markers for each trait was determined as follows (Utz et al. 2000): where R 2 adj is the adjusted R 2 and H 2 is the broad-sense heritability of the trait. The adjusted R 2 was estimated by fitting a multiple linear regression model as follows: where y is the BLUEs and snp i , snp j , and snp k are the SNP markers identified by ST-GWAS with threshold of snp i > snp j > … > snp k (Gaikpa et al. 2020). Individual p G values were calculated using the sum of squares (SS) from the analysis of variance of the linear model (Würschum et al. 2015) as: where SS snp is the sum of squares of individual SNP and SS total is the total sum of squares of the model. The sum of squares was estimated using the SS type II to avoid changes in the estimates due to SNPs order in the model. Additive effects were estimated by fitting individual SNP in the linear model one at a time as described by Gaikpa et al. (2020).
Additionally, pleiotropic loci between FHB severity, PH and AR were investigated by performing a multi-trait GWAS (MT-GWAS) using the multi-trait mixed-model (MTMM) method proposed by Korte et al. (2012). MT-GWAS was implemented in ASRemL-R package v.4.1.0 between FHB severity and PH (FHB vs PH), FHB severity and AR (FHB vs AR) and PH and AR (PH vs AR). For each pair of traits, the MTMM was: where y 1 and y 2 are the phenotypes of the first and second traits, respectively; S i is a vector of 1 for all values belonging to the ith trait; µ i is a vector of fixed effects including the mean of the ith trait; X is the vector of SNP markers; β is the vector of common effects between the two traits; α is the vector of interaction effects; and υ is a random variable capturing both the error and genetic random effects such as the kinship matrix (K). To identify causal loci with common effects (COM) as well as opposite/interaction effects (IE) on each pair of traits, we performed three generalized least squares (GLS) F tests (Korte et al. 2012). Firstly, the full model was tested against a null model where β = 0 and α = 0 in order to identify the combination of both COM and IE loci. Secondly, COM loci were detected by testing the model with α = 0 against the null model. Thirdly, the model with β = 0 was tested against the full model to identify IE loci. Significant effects were identified based on the same exploratory and Bonferroni-corrected threshold used in ST-GWAS. This helped to further reduce potential false positives due to inflation of p values from F tests as found by Merrick et al. (2022). The genotypic correlation between traits pertaining to each pleiotropic locus (r snp ) was estimated by fitting the significant markers as fixed effects in the bivariate model described in Eq. (5). Markers were fitted consecutively one after the other starting from the highest p value to the lowest; that is, threshold of snp i > snp j > … > snp k . In that order, the genetic correlation attributable to the jth SNP was the difference between the correlation of the model after fitting the ith SNP and the one detected after the jth SNP. Furthermore, the proportion of genetic correlation explained by each pleiotropic locus was calculated as: where p Ci is the proportion of genotypic correlation explained by the ith SNP and r g is the total genotypic correlation obtained from Eq. (5). For both ST-GWAS and MT-GWAS, the genomic relationship matrix was estimated using the kinship algorithm, VanRaden (2008) in GAPIT v.3.1.0. The kinship matrix (K) was fitted as covariate variable in each GWAS model. Quantile-quantile (Q-Q) plots were used to evaluate the power of the model in controlling false positives and negatives. Q-Q plots having a straight line, close to 1:1 on the diagonal, with a sharp upward deviated tail, indicate 1 3 a good control of both false positives and negatives by the model and confirm the existence of true MTAs (Kaler et al. 2020). Q-Q and Manhattan plots were drawn using the R package CMplot v.4.0.0 (Yin et al. 2021). Moreover, linkage disequilibrium (LD) patterns were also investigated by calculating R 2 values between pairs of significant markers to identify chromosomal positions of unmapped SNP markers. Markers were considered to be in LD when R 2 ≥ 0.70.

Genomic prediction and marker-assisted selection
Genomic prediction analysis was performed using the mixed linear ridge regression model in the rrBLUP package v.4.6.1 (Endelman 2011). The model (GP 1 ) included all markers as random effects and was fitted as follows: where Y is the vector of BLUEs; β is a vector of fixed effect such as the grand mean; Z is the design matrix; u ~ N(0, A 2 ) is the vector of random markers effects; A is a relationship matrix and the residuals are normal with constant variance; and ε is a vector of residual errors. GP 1 model was validated using the fivefold cross-validation approach where the training and validation sets were constituted by splitting the phenotypic and genotypic data into five sets, consisting of 80-81 genotypes each (Gaire et al. 2022). Sampling was repeated 50 times, and GP 1 model was run for each set to determine the genomic prediction ability (r MG ) and accuracy (r MG/H ) for each of the traits. Genomic prediction ability was the Pearson correlation coefficient between predicted and observed phenotypes for each of the five validation sets (Ould Estaghvirou et al. 2013). Moreover, the prediction accuracy was determined for each validation set as described by Ould Estaghvirou et al. (2013): where H 2 is the broad-sense heritability of the trait. rMG and rMG/H values were averaged over the five validation sets to obtain the prediction ability and accuracy of each GP 1 model. Additionally, we evaluated the potential of markerassisted selection (MAS) using only significant markers from ST-GWAS which explained at least 5% of the genetic variation. MAS prediction was done based on the markers effects deducted from the GP 1 model, and the prediction accuracy was calculated.

Estimation of genomic background effect for PH, FHB severity and AR and selection of resistant genotypes with Rht-D1b
Based on the ST-GWAS, we evaluated the contribution of GB to FHB severity, PH and AR using the additive effects and estimates of total genetic variation (p G values) explained by GB and Rht markers. p G values were estimated separately for GB and Rht markers by fitting a linear regression model and using the adjusted R 2 as described in Eq. (15). In addition, a second genomic prediction model (GP 2 ) was built and the prediction ability (r MG ) and accuracy (r MG/H ) were evaluated following the fivefold approach. Training and validation sets were the same as described for GP 1 . In GP 2 , significant MTAs for PH on chromosomes 4D and 6A were used as covariates, and the prediction accuracy was estimated for each trait based on small-effect markers only. GP 2 was fitted as follows: where Y is the vector of BLUEs; β is the vector of fixed effects of the covariates; X and Z are the design matrices; u ~ N(0, A 2 ) is the vector of random markers effects; A is a relationship matrix and the residuals are normal with constant variance; and ε is a vector of residual errors. Furthermore, we selected the ten best genotypes with low FHB severity and Rht-D1b allele based on the calculation of genetic estimated breeding values (GEBV) using the effects of random markers from the GP 2 model as: where Y 0 and Z 0 are the vectors of GEBV and design matrix of the genotypes, respectively. The correlation between the GEBV-based approach and stacking of additive effects (SAE) of GB markers from ST-GWAS was also determined in order to check the relative effectiveness of GWAS for exploiting genomic background. SAE was calculated only for GB markers which explained at least 5% of genetic variation.

Considerable genetic variation was found for FHB severity and morphological traits
Genetic coefficients of variation (CV G ) were low for DTF, HS, EL, NS and ED and moderate to high for PH, FHB severity and AR (Table 1). The higher CV G observed for AR was due to the considerable existing genetic variation among the genotypes. The coefficient of variation due to was low for all traits (Table 1). Genetic variances ranged from 0.03 to 6.9 for DTF, HS, EL, NS and ED and were higher (97.4-489.4) for PH, FHB and AR. The genotype-environment interaction (GEI) was significant and generally lower than the genetic variance. Depending on the trait, GEI variance was 1.5-22.2-fold lower than the genetic variance. Broad-sense heritability estimates were very high (> 0.9) for all traits.

Correlations and path coefficients depict the contribution of morphological traits to FHB severity
Density plots of BLUEs across environments showed nearly normal distribution for all traits (Fig. 2). Linear relationships with significant phenotypic correlations were observed among traits. FHB severity was positively and highly correlated with AR (r p = 0.70) and negatively and moderately correlated with PH (r p = − 0.62). The outlying point in the scatter plot FHB severity/AR belongs to the cleistogamic cultivar "Anapolis" that has a very high AR (99.6%) and a low FHB severity (20.9%). PH was negatively and moderately correlated with AR (r p = − 0.53), while NS showed positive and low correlation with EL (r p = 0.39) and ED (r p = 0.37). Highly negative phenotypic correlations were detected between EL and ED (r p = − 0.71) (Fig. 2). Phenotypic correlations between FHB severity and EL and NS were significant, but low. At the genetic level, the correlations were negative and moderate between FHB severity and PH (r g = − 0.64), PH and AR (r g = − 0.55) and negatively high between ED and EL (r g = − 0.70). Genotypic correlations were positively high between FHB severity and AR (r g = 0.74) and low between NS with EL (r p = 0.38) and ED (r p = 0.37) ( Table 2). The causal system formed by PH, AR, EL and NS explained 67.04% of the variation of FHB severity (Table 3). Path effect due to residual factors was 0.56. The direct effect of AR (0.57) was 3.3-fold higher than its indirect paths and similar to the residual effect. The indirect effect of PH via AR (− 0.31) was slightly higher than its direct effect (− 0.28), but both effects were lower than the residual effect. Direct and indirect effects of EL were considerably lower than the residual effect, while effects of NS were not significant. Inflation due to multicollinearity among PH, AR, EL and NS was negligible (VIF < 2), showing that estimates of direct and indirect effects and coefficient of determination were accurate (Table 3).

Several marker-trait associations (MTAs) were detected by genome-wide association (GWA) studies
The first two principal coordinates (PCs) indicated that the 401 genotypes included in the study were not structured (Fig. S1). Therefore, we corrected in the GWAS for relatedness among genotypes using the kinship matrix and did not include PCs. The ST-GWAS was performed to depict the genetic architecture of FHB severity, PH and AR (Table S2). Considering all traits, a total of 25 significant MTAs were identified, of which seven were specific to ST-GWAS 1 , seven to ST-GWAS 2 and eleven common to both models (Fig. 3a, b). The use of Rht markers as covariate variables in ST-GWAS 2 identified additional small-effect MTAs for all traits. Particularly, the number of small-effect MTAs increased from four to seven for AR, while two new MTAs were identified for FHB severity. In total, eight, eleven and nine MTAs were identified for PH, FHB severity and AR, respectively (Table S2). One MTA was common to PH and FHB severity, one to PH and AR and two to FHB severity and AR. MTAs explained 69.7, 60.7 and 44.2% of total genetic variation for PH, FHB severity and AR, respectively. Markers rs20873 on chromosome 4D linked to Rht-D1 and rs10110 on chromosome 6A linked to Rht24 were major MTAs for PH (Fig. 3a, b), which explained 45.2 and 10.8% of genetic variation, respectively. In addition, rs20873 conveyed the largest reduction effect of − 7.02 cm on PH. For FHB severity, major MTAs were rs20873, rs5192 on chromosome 7B and rs3647 on chromosome 5A. Likewise, rs20873 was the major MTA identified for AR, explaining 22.7% of genetic variation. rs20873 increased FHB severity and AR by about 6 and 12%, respectively. The false discovery rates were very low (FDR < 0.05) for identified MTAs. The MT-GWAS identified a total of six pleiotropic loci distributed on chromosomes 2A, 4D, 5A, 6B and 7B. All six were detected between FHB severity and PH (Table 4a), four between FHB and AR (Table 4b) and three between PH and AR (Table 4c). In comparison with the ST-GWAS from BLINK, the marginal trait from the MT-GWAS identified three additional MTAs on chromosomes 4D, 6B and 5A for FHB severity, one on chromosome 5A for PH and one on chromosome 4D for AR. However, the MT-GWAS identified three MTAs for FHB severity, one for PH and three for AR, which were also detected by ST-GWAS. Three out of the six MTAs detected for FHB versus PH exhibited positive contribution to the genotypic correlation and thus exerted common effects on FHB severity and PH (Table 4a). Inversely, the remaining three MTAs showed negative genotypic correlations, revealing their opposite effects on both traits. Marker rs9660 on chromosome 4D was associated with FHB severity only, but it exerted an opposite effect on the two traits. Similarly, all common effects MTAs were associated with FHB severity only. Common effect MTAs explained 6.3-26.6% of the genotypic correlation between FHB severity and PH, while 3.1-6.3% of the correlation was attributable to interaction effects MTAs (Table 4a). For FHB versus AR, all four MTAs identified by the full model had common effects with 1.4-14.9% of the genotypic correlation explained between FHB severity and AR (Table 4b).
Although marker rs5192 on chromosome 7B was specific to FHB severity, it had a common effect on both traits. For PH versus AR, two out of the three MTAs identified by the full model had opposite effects with 5.5-56.4% of the genotypic correlation explained between the two traits (Table 4c). Marker rs7788 that was specific to AR had a common effect on both PH and AR with 10.9% of the genotypic correlation explained.
In both ST-GWAS and MT-GWAS, the inspection of Q-Q plots of the two GWAS revealed a good control of false positives and negatives, demonstrating the existence of true MTAs controlling the genetic architecture of each trait as well as their interactions (Figs. S2-S5). The linkage disequilibrium analysis revealed strong linkage (R 2 = 0.9) between markers rs19377 and rs13973 located on chromosome 5B for PH (Fig. S6). Very low R 2 values were observed among all other MTAs.

Significant differences were observed among groups of genotypes
NoRht, Rht24b, Rht-D1b and Rht24b + Rht-D1b represented, respectively, 23.7, 27.7, 14.9 and 33.7% of our materials (Table 5). Furthermore, considerable genetic variation was observed within each group for PH, FHB severity and AR (Table 5). For PH, the genetic variance was 2.5-6.1-fold higher in NoRht compared with other groups. The lowest genetic variance was observed in Rht-D1b for PH. However, in FHB severity, the genetic variance was 1.3-1.3-fold higher in Rht-D1b than in NoRht and Rht24b. The genetic variance was also 1.5-1.6-fold higher in Rht24b + Rht-D1b compared with NoRht and Rht24b which had similar variation for FHB severity. Genetic variance was high and similar among groups for AR. Genotype-environment interactions were significant and relatively low for the three traits in all groups. Comparative analysis among the four groups revealed statistically significant differences for all traits (Fig. 4). Genotypes were on average shorter in Rht24b + Rht-D1b and taller in NoRht. In Rht24b, genotypes were taller than those in Rht-D1b and Rht24b + Rht-D1b. Contrary to PH, FHB severity was on average lower in NoRht and Rht24b than in Rht-D1b and Rht24b + Rht-D1b. FHB severity was statistically identical between NoRht and Rht24b, and Rht-D1b and Rht24b + Rht-D1b, respectively. Similar to FHB severity, the average AR was lower in NoRht and Rht24b compared with Rht-D1b and Rht24b + Rht-D1b, respectively. No significant difference was found for AR between NoRht and Rht24b, as well as Rht-D1b and Rht24b + Rht-D1b. The reduction effect of Rht alleles on PH was moderately correlated with FHB severity and AR. Highly negative reduction effect was associated with high FHB severity and AR.

Genomic prediction (GP) exhibited high prediction ability (rMG) and accuracy (rMG/H) for FHB severity, PH and AR
Genomic prediction ability and accuracy averaged over the five validation sets were greater than 0.7 depending on the trait (Table S3). FHB severity and AR showed similar and lower prediction accuracy than PH. Genotypes with Rht-D1b and Rht24b alleles were distributed across all the validation sets, showing a good representation of the genetic variation in each set (Table S4). This guaranteed the accurate estimation of both r MG and r MG/H , resulting in low standard errors of the estimates across validation sets (Table S3). In contrary, the prediction accuracy of marker-assisted selection based on SNP markers with p G ≥ 5% was moderate for all traits, with the lowest value (0.4) observed in AR (Fig. 5). GP 1 including all available markers increased the prediction accuracy by about 0.2, 0.3 and 0.2 for FHB severity, AR and PH, respectively.

Genomic background affects PH, FHB severity and AR
The contribution of GB to the genetic variation was different among the three traits (Fig. 6). The total proportion of genetic variation explained by GB markers was higher for FHB severity than PH and AR (Fig. 6a). Contrary to PH, GB markers had higher contribution to the genetic variation of FHB severity and AR than Rht markers (Fig. 6a). Ten, eight and six GB markers were found for FHB severity, AR and PH, respectively (Fig. 6b). Non-Rht markers which explained at least 5% of the genetic variation were rs5192, rs3647 and rs13165 for FHB severity, rs3629 and rs8776 for AR and rs13233 for PH (Fig. 6b). Reduction effects exerted by GB markers were relatively high for all traits (Fig. 7a-c). The cumulative reduction effect (CRE) calculated from markers with p G ≥ 5% was − 5.7 cm, − 11.2% and − 14.3% for PH (Fig. 7a), FHB severity (Fig. 7b) and AR (Fig. 7c), respectively. CRE was about twofold lower than the additive effects of Rht markers for PH. Inversely, the CRE exerted by GB markers on FHB severity and AR was 1.3-1.9-fold higher than the increase effect of Rht markers (Fig. 7b, c). Results of the second genomic prediction (GP 2 ) implemented using Rht markers as covariates showed moderate prediction accuracy for the three traits (Table S5). This showed that the estimated contribution of GB was higher in the three traits with GP 2 compared with GWAS. The highest accuracy was observed for AR, while the lowest value was obtained in PH, confirming the trend revealed by the GWAS. Standard errors were also low for GP 2 , indicating a high consistency of the prediction accuracy among the five validation sets.

Genotypes with low FHB severity in the presence of Rht-D1b were selected using genomic background
Depending on Rht groups, GEBV from GP was positively correlated (r = 0.64-0.77) with stacking of additive effects (SAE) of GB markers with p G ≥ 5% (Fig. 8). However, several genotypes with same SAE (i.e. − 14.42) exhibited different GEBV. The phenotypic performances of the ten best genotypes selected based on GEBV from the genomic prediction (GP 2 ) in Rht-D1b and Rht24b + Rht-D1b groups are summarized in Table 6. The negativity and size of the effects of GB markers as measured by GEBV describe the degree of resistance GB for FHB severity. Thus, highly negative GEBV for FHB severity indicates higher resistance GB, whereas positive GEBV indicates low GB or susceptibility GB. Genotypes with the lowest observed FHB severity exhibited negatively high GEBV for FHB severity carrying Rht-D1b allele only (e.g. Faktor and Anapolis), or the combination Rht24b + Rht-D1b alleles (e.g. Kranich and Kamerad) ( Table 6). More resistant genotypes within each Rht group also showed negatively high GEBV with lower observed AR, except Anapolis which had a highly positive GEBV for AR and the highest observed AR at all. Taller genotypes had positive GEBV for PH.

Discussion
Through a combined correlation and path analysis, we provided a first insight into the nature and magnitude of the complex interactions between FHB resistance and morphological traits including plant height and anther retention. ST-GWAS, MT-GWAS and genomic prediction were implemented to depict the genetic architecture of FHB severity, anther retention and plant height and understand the effect of alternative Rht genes on FHB resistance. Additionally, multi-trait GWAS was conducted to identify positively and negatively pleiotropic loci controlling complex interactions among traits. The ST-GWAS and genomic prediction helped to evaluate for the first time the contribution of genomic background to each trait and its potential to improve FHB resistance in wheat genotypes with semi-dwarfing allele Rht-D1b.

Existence of high path effects set criteria for effective indirect selection for FHB resistance
FHB severity exhibited high positive correlation with anther retention and moderate negative correlation with plant height at both phenotypic and genotypic levels. Similar observations were reported by Buerstmayr and Buerstmayr (2015); Steiner et al. (2019b) and Ruan et al. (2020). Buerstmayr and Buerstmayr (2015) reported that FHB severity was positively correlated with anther retention (0.63) and negatively with plant height (− 0.39). Previous investigations on anther extrusion, the opposite of anther retention, also highlighted negative correlation of − 0.45 to − 0.64 (Lu et al. 2013) and − 0.55 to − 0.74 (Xu et al. 2020). Recently, Nannuru et al. (2022) also reported negative moderate correlations between FHB severity with anther extrusion (− 0.48) and plant height (− 0.43) in Nordic spring wheat. Our results confirmed that anther retention was negatively correlated with plant height, indicating the existence of shared genetic control between the two traits in winter wheat in contrary to Nannuru et al. (2022) who found a lower correlation (0.16) between plant height and anther extrusion. The differences in the strength of the correlations across studies are an evidence that interactions among the traits are also affected by environment and breeding material characterized by diverse genetic factors.
The causal system formed by morphological traits including anther retention and plant height explained up to 67% of the phenotypic variation of FHB severity in this study. This confirms the importance of morphological traits in the passive resistance mechanism against FHB disease in winter wheat (Mesterházy 1995;Buerstmayr et al. 2020;Xu et al. 2020). The effect due to residual factors was high (0.56) and can be attributed to specific unshared local genetic factors which had considerable influence on FHB severity. The existence of high direct path of 0.57 between FHB severity and anther retention indicates that any increase by 1% in the standard deviation of anther retention implies a direct increase of 0.57% in the standard deviation of FHB severity independently from other traits. This finding exhibits anther retention as a major indicator trait to be included in multipletrait breeding strategies to aim for FHB-resistant cultivar development in wheat. Particularly, when the breeder's interest is not in other traits, an indirect selection using anther retention is likely to yield FHB-resistant cultivars (Fernandes et al. 2018;Moreno-Amores et al. 2020). However, path effects between plant height and FHB severity were lower than the residual and the direct effect of anther retention, indicating that plant height as a sole trait does not have considerable contribution to the variation of FHB severity. It is worth noting that the indirect path effect of plant height via anther retention was higher than its direct effect on FHB severity, and this can be explained by the association between plant height and anther retention and underlying genetic factors. A decrease of 1 cm in the standard deviation of plant height in a cultivar would cause a direct increase of 0.28% in the standard deviation of FHB severity   and an indirect increase 0.31% via anther retention. This demonstrates firstly that the moderate correlation detected between plant height and FHB severity was mainly due to anther retention, and secondly that plant height per se has a high impact on FHB severity, even when the experiments are inoculated from above like in this study. As implication, plant height should therefore be considered simultaneously with anther retention in a multiple-trait selection approach for the development of semi-dwarf FHB-resistant cultivars.

Combined ST-GWAS and MT-GWAS revealed two major Rht genes and novel pleiotropic loci for plant height, FHB severity and anther retention
Frequency distributions (Fig. 2) already indicated that all traits analysed here were quantitatively inherited, i.e. by many loci of varying effects. Indeed, a total of eight loci that control plant height were identified, with corresponding  (Würschum et al. 2015(Würschum et al. , 2017Herter et al. 2018;Khadka et al. 2021;Tian et al. 2022), respectively, explaining 45.31% and 10.81% of the total genetic variation of plant height. This supports findings of Würschum et al. (2017), Herter et al. (2018 and Tian et al. (2017), who found that Rht24 was the second most important Rht gene in Central European and worldwide commercial wheat breeding programmes. Rht-D1 and Rht24 are, respectively, gibberellin-insensitive (De Velde et al. 2021) and gibberellin-sensitive (Tian et al. 2022) genes responsible for reduced plant height in wheat. The other genomic regions were also enriched with several medium-and small-effect QTL which contribute to plant height. Those small-effect QTL represent the genomic background which explains existing genetic variation of plant height within Rht groups. In addition to Rht-D1, a larger number of small-loci controlling FHB severity were identified on other genomic regions such as 2A, 2B, 4A, 5A, 6A, 6B and 7B, which were considered as genomic background contributing the variation of FHB severity among semi-dwarf genotypes. Several FHB resistance QTL were also reported on chromosomes 2A (He et al. 2016a;Gadaleta et al. 2019), 2B (Sari et al. 2018;Ollier et al. 2020), 4A (He et al. 2016a;Zhang et al. 2021;Nannuru et al. 2022), 5A (Steiner et al. 2019a;Ruan Fig. 8 Scatter plot showing the strength of relationship between stacking of additive effects (SAE) of genomic background (GB) markers from single-trait genome-wide association study (ST-GWAS) and genomic estimated breeding values (GEBV) from genomic prediction (GP). SAE was estimated based on GB markers with p G ≥ 5%. ***significant at p < 0.001  Sari et al. 2020;Nannuru et al. 2022), 6A , 6B (Cuthbert et al. 2007;Ollier et al. 2020) and 7B (Sari et al. 2018;Ollier et al. 2020). In addition, FHB resistance QTL on chromosomes 5A and 7B explained more than 7% of the genetic variation each and can be considered as important locally adapted QTL to improve FHB resistance in European winter wheat. The QTL on chromosome 5A was linked to marker (AX-158558712) at position 521 Mbp (Table S6). This falls within the QTL region Qfhb. nmbu.  which was reported through a meta-QTL analysis by Venske et al. (2019) and recently detected in the European wheat panel by Nannuru et al. (2022). For anther retention, in addition to Rht-D1b, several medium-and small-effect QTL were detected in genomic regions including chromosomes 2A, 3A, 4A, 5A and 7A. Previous studies also identified QTL for anther retention on chromosomes 3A (Muqaddasi et al. 2019), 4A (Buerstmayr and Buerstmayr 2015) and 5A (Buerstmayr and Buerstmayr 2015;Steiner et al. 2019a;Sari et al. 2020;Xu et al. 2020). These QTL, representing the genomic background, contribute to the adjustments of anther retention, particularly within Rht genes groups.
In contrary to Rht24b, Rht-D1b exerted an opposite or interaction effect on both plant height and FHB severity. Rht-D1b explained more than 26% of the observed negative genotypic correlation between the two traits. The linkage disequilibrium analysis revealed an absence of close association between QTL identified for the different traits. This exhibits Rht-D1 as a major gene conveying a negatively pleiotropic effect on FHB resistance. As indicated by Raherison et al. (2020), alleles of negatively pleiotropic loci exert favourable effects on one trait and unfavourable effects on the other trait depending on the breeding objectives. In commercial wheat breeding, the reduction effect of Rht-D1b on plant height represents a favourable effect, while its increase effect on FHB severity or susceptibility is perceived as an unfavourable effect. Several studies also found that Rht-D1b was highly associated with FHB susceptibility (Srinivasachary et al. 2008;Mao et al. 2010;Lu et al. 2013;Buerstmayr and Buerstmayr 2016;He et al. 2016b;Prat et al. 2017;Hales et al. 2020;Zhang et al. 2020). Two other interaction effect loci were also identified in other regions of chromosomes 4D and 5A between plant height and FHB severity. The second causal locus on chromosome 4D with linked marker Bob-White_rep_c48828_217 at position 195 Mbp was associated with FHB severity only, demonstrating the existence of mediated negatively pleiotropic loci that contribute to the negative correlation between FHB resistance and plant height. Mediated pleiotropy as described by Hackinger and Zeggini (2017) represents a situation where a causal locus controls one trait, which in turn causes a second trait, resulting in significant association between the two traits as detected by correlation/covariance analysis. Moreover, our results also revealed that the small-effect loci identified for FHB severity on chromosomes 2A and 7B by the ST-GWAS and another locus on chromosome 6B conveyed a common effect on both FHB resistance and plant height. Although associated with FHB severity only, these QTL exhibited positive contribution to the interactions between the traits with about 16% of the genotypic correlation explained. Ghimire et al. (2022) also reported five stable and pleiotropic QTL associated with FHB resistance and related traits such as DON accumulation and plant height in red winter wheat. Similarly, Schulthess et al. (2017) reported the existence of genomic regions inducing pleiotropic effects on grain yield and related traits in wheat. Our study revealed the existence of positively pleiotropic loci which changes (increase or decrease) FHB severity and plant height in the same direction. These positively pleiotropic QTL can be exploited in breeding programmes to reduce FHB severity (i.e. improve FHB resistance) and plant height simultaneously. Most importantly, the QTL on chromosome 7B linked to marker AX-158601365 at position 661 Mbp, explaining more than 7% of the genetic variation of FHB severity and exerting a positively pleiotropic effect on both FHB severity and anther retention represents a novel promising QTL which could be integrated into multiple-trait breeding strategy for higher FHB resistance in European winter wheat. Breeder-friendly KASP markers can be developed to facilitate the integration of this QTL into marker-assisted selection.
The high positive correlation and direct path effect observed between FHB severity and anther retention were controlled by the pleiotropic loci identified between FHB severity and plant height, particularly Rht-D1 and the three QTL distributed on chromosomes 4D, 2A and 7B. All pleiotropic loci exerted common effects on FHB resistance and anther retention with about 24% of genotypic correlation explained in total. This firstly supports the existence of shared small-effect QTL between FHB resistance and anther retention (Lu et al. 2013) and secondly demonstrates that shared QTL may have positively pleiotropic effects on the two traits. Moreover, this positive pleiotropy can be either direct or mediated as the causal pleiotropic locus on chromosome 7B was specific to FHB severity, while others were associated with both traits. Existence of positive pleiotropy between FHB severity and anther retention offers the opportunity of utilizing genomics-assisted breeding to improve both FHB resistance and anther retention in winter wheat. Considering that the identified pleiotropic loci were smalleffect QTL, except Rht-D1, the efficient exploitation of positive pleiotropy for improving FHB severity and anther retention could be achieved by implementing multi-trait genomic prediction as reported by Gaire et al. (2022) for FHB-related deoxynivalenol accumulation in wheat. Furthermore, the moderate negative correlation observed between plant height and anther retention was controlled by Rht-D1 and the two pleiotropic loci on 2A and 5A. Rht-D1 and QTL on 5A exhibited negatively pleiotropic effects on plant height and anther retention, while the QTL on 2A exerted a mediated positively pleiotropic effect on the two traits. Rht-D1 explained more than 56% of observed correlation between these traits, indicating that Rht genes had a major impact on anther retention and supporting that Rht-D1b was associated with high anther retention as recently reported by Buerstmayr and Buerstmayr (2022).

Rht24b has no effect on FHB severity and anther retention
The combination of Rht-D1b and Rht24b was more frequent in our cultivars as previously reported within European wheat germplasm (Würschum et al. 2017;Tian et al. 2019). The similar FHB severity observed between Rht24b and genotypes with tall alleles (NoRht) is an indication that Rht24b did not affect FHB resistance (Herter et al. 2018;Miedaner et al. 2022). Similarly, our results also revealed that Rht24b did not significantly contribute to anther retention. Accordingly, no marker was found that co-segregates with the other morphological traits (Fig. 3). This offers the possibility of developing semi-dwarf genotypes with improved FHB resistance and low anther retention. However, the effect of Rht24b allele (− 3.42 cm) on plant height was about twofold lower than the effect of Rht-D1b (− 7.85 cm, Fig. 7a), showing that Rht24b exerted a lower reduction effect on plant height compared with other Rht genes (Tian et al. 2017;Herter et al. 2018). In addition, the high genetic variation for FHB severity observed among genotypes with Rht24b indicates the random existence of FHB-susceptible genotypes within this group. With this, alternative sources of resistance including locally adapted loci could also be explored to develop FHB-resistant genotypes within each Rht gene group, depending on breeding objectives.

Genomic background has the potential to improve FHB resistance in genotypes with Rht-D1b
Considerable genetic variation was observed for FHB severity and anther retention among genotypes with Rht-D1b and Rht24b + Rht-D1b alleles, demonstrating that there is room for the development of FHB-resistant cultivars with the Rht-D1b allele. This large genetic variation in semi-dwarf genotypes can be attributed to the effects of genomic background distributed across several genomic regions. A key implication is that the genomic background has the ability to counterbalance the negative effect of Rht-D1b on FHB resistance and anther retention as explained by Buerstmayr and Buerstmayr (2022) who demonstrated by backcrosses with four near-isogenic lines (NILs) that the background resistance of the lines reduced efficiently the effect of semidwarfing alleles on FHB severity in spring wheat. Brar et al. (2019) also found that the genomic background and epistatic interactions had a significant impact on the expression of FHB resistance in hard red spring wheat. The more negative the genomic background effect, referring to higher resistance genomic background, the lower the FHB severity and anther retention. This offers a great opportunity to exploit the genomic background for improving FHB resistance in dwarf genotypes.
Effects of the genomic background on FHB severity and anther retention can be efficiently evaluated by calculating the genomic estimated breeding values (GEBV) of genotypes based on the effects of all loci, except Rht genes in our case (Bonnett et al. 2022). Several studies have suggested the stacking of additive effects or favourable alleles of QTL from genome-wide association study as an efficient approach to improve FHB resistance and related traits (Miedaner et al. 2006;Sidhu et al. 2020;Ghimire et al. 2022;Nannuru et al. 2022). However, our results showed that despite the existence of moderate to high correlations (r > 0.6) between the GEBV-based and stacking of additive effects (SAE) approaches, many genotypes with similar additive effects from GWAS exhibited different GEBV and FHB severity. This shows the limitations of SAE to accurately discriminate among genotypes based on their genomic background contrary to the GEBV-based approach and hence confirms the superiority and efficiency of the genomic prediction over fixed effects selection methods as reported by several previous studies (Juliana et al. 2017(Juliana et al. , 2022Sandhu et al. 2021). In addition, the higher accuracy observed for the genomic prediction compared with marker-assisted selection based on the effects of QTL with p G ≥ 5% from the ST-GWAS (Fig. 5) is another indication that genomic prediction would help to better exploit genomic background than GWAS. This poor performance of GWAS can be explained by the quantitative nature of FHB resistance and anther retention Ruan et al. 2020), and that genomic background is constituted of several medium-and smalleffect QTL. However, the genome-wide association study (GWAS) can identify only a small number of those minor QTL, not regarding the remaining which are also part of the genomic background and could have together a significant contribution to the phenotype. Moreover, using the GEBVbased approach, the best two genotypes with Rht-D1b had a FHB severity of 20.3 and 20.8%, respectively ( Table 6). The best two genotypes without dwarfing alleles (Carimulti, Helmond, Table S1) had a FHB severity of 12.4 and 12.8%, respectively. Hence, Rht-D1b still has a penalty for FHB resistance of about 40% of increase in disease severity compared to genotypes without this semi-dwarfing allele. This illustrates the necessity to use alternative Rht alleles including Rht24b, boosted by a short genomic background or other FHB-neutral Rht dwarfing alleles.
Furthermore, the GEBV approach would select Anapolis as one of the best genotypes which exhibited the highest anther retention (99.6%). This is a typical cleistogamous genotype which possesses a structural floral barrier for pathogens, resulting in low FHB severity (Kubo et al. 2010;Tang et al. 2020;Buerstmayr et al. 2021;Zajączkowska et al. 2021). Similar observations were made in barley by Kawada and Kubo (2008) who found that cleistogamous genotypes had high FHB resistance and low mycotoxin accumulation.

Concluding remarks
A clear understanding of the complex interactions between FHB severity and morphological traits and their genetic architecture can significantly contribute to addressing the needs for semi-dwarf wheat genotypes with high FHB resistance. The existence of high direct and indirect path effects between FHB severity and morphological traits demonstrates that multiple indirect trait selection for FHB severity has a great potential and should always integrate anther retention and plant height as important secondary traits. Positively direct and/or mediated pleiotropic loci controlling complex interactions between traits could be integrated into breeding programmes for an efficient multipletrait selection for higher FHB resistance in wheat. Genomic background (GB) explained a high proportion of the genetic variance of FHB severity and anther retention and has the potential to increase FHB resistance in genotypes with Rht-D1b. Depending on the breeding goals, the development of FHB-resistant cultivars should consider Rht24b which was not linked to FHB susceptibility and exploit the GB for higher FHB resistance to counterbalance the negative effect of Rht-D1b. Similarly, PH could be further reduced in the Rht24 group by the exploitation of GB. Strategies to efficiently exploit existing interactions among traits and the GB in breeding programmes to develop FHB-resistant cultivars with Rht-D1b should include the selection of: (i) taller genotypes within the Rht-D1b group as long as they have a high (negative) FHB-related resistance GB (e.g. Faktor, Mikon, Table 6), (ii) genotypes having the lowest anther retention, exploiting the high direct path and positively pleiotropic loci between FHB severity and anther retention, and/or (iii) genotypes with anther retention higher than 95% indirectly taking advantage of the contribution of cleistogamy to passive FHB resistance (e.g. Anapolis). Such genotypes could be most efficiently detected by using genomic estimated breeding values for FHB severity.
Human and animal rights This research does not contain any studies on humans or animals.
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/.