Genetic mapping and genomic prediction of sclerotinia stem rot resistance to rapeseed/canola (Brassica napus L.) at seedling stage

GWAS detected ninety-eight significant SNPs associated with Sclerotinia sclerotiorum resistance. Six statistical models resulted in medium to high predictive ability, depending on trait, indicating potential of genomic prediction for disease resistance breeding. The lack of complete host resistance and a complex resistance inheritance nature between rapeseed/canola and Sclerotinia sclerotiorum often limits the development of functional molecular markers that enable breeding for sclerotinia stem rot (SSR) resistance. However, genomics-assisted selection has the potential to accelerate the breeding for SSR resistance. Therefore, genome-wide association (GWA) mapping and genomic prediction (GP) were performed using a diverse panel of 337 rapeseed/canola genotypes. Three-week-old seedlings were screened using the petiole inoculation technique (PIT). Days to wilt (DW) up to 2 weeks and lesion phenotypes (LP) at 3, 4, and 7 days post-inoculation (dpi) were recorded. A strong correlation (r = − 0.90) between DW and LP_4dpi implied that a single time point scoring at four days could be used as a proxy trait. GWA analyses using single-locus (SL) and multi-locus (ML) models identified a total of 41, and 208 significantly associated SNPs, respectively. Out of these, ninety-eight SNPs were identified by a combination of the SL model and any of the ML models, at least two ML models, or two traits. These SNPs explained 1.25–12.22% of the phenotypic variance and considered as significant, could be associated with SSR resistance. Eighty-three candidate genes with a function in disease resistance were associated with the significant SNPs. Six GP models resulted in moderate to high (0.42–0.67) predictive ability depending on SSR resistance traits. The resistant genotypes and significant SNPs will serve as valuable resources for future SSR resistance breeding. Our results also highlight the potential of genomic selection to improve rapeseed/canola breeding for SSR resistance.


Introduction
Rapeseed/canola (Brassica napus L., genomes = AACC, 2n = 4x = 38) is an amphidiploid Brassica species and the second largest cultivated oilseed crop in the world after soybean (USDA Foreign Agricultural Service 2021). Sclerotinia stem rot (SSR), caused by the necrotrophic plant pathogenic fungus Sclerotinia sclerotiorum (Lib) de Bary, is one of the most economically important diseases affecting rapeseed/ canola that significantly limits worldwide rapeseed/canola production (Boland and Hall 1994;Bolton et al. 2006). The yield losses due to this pathogen vary from 10-80% from year to year depending on the disease development environments (Del Río et al. 2007;Wu et al. 2016). However, in the USA, each unit increase in SSR incidence imposes 0.5-0.7% loss in canola seed yields (Del Río et al. 2007;Koch et al. 2007). Moreover, SSR affected plants often tend to have reduced oil content, and inferior oil quality due to the changing of oil's fatty acid profile (McCartney et al. 1999;Sharma et al. 2015).
To manage the associated risk of SSR disease, growers primarily depend on the use of conventional rotation with non-host crop species and chemical controls which are neither completely effective nor economically and environmentally feasible (Derbyshire and Denton-Giles 2016;Roy et al. Communicated by Jacqueline Batley. 2021). Therefore, breeding for durable SSR resistant varieties would be a more economically feasible, environmentfriendly, and sustainable strategy to manage this disease. However, no accessions conferring high level of resistance or complete immunity to S. sclerotiorum have been identified over the last three decades of investigation (Zhao et al. 2004;Bradley et al. 2006;Yin et al. 2010). Thus, the current breeding strategy for improved SSR resistance is solely dependent on the utilization of partially resistant germplasm. Therefore, it is crucial to screen a worldwide collection of diverse genotypes with an appropriate screening method to identify genetically resistant genotypes to improve SSR resistance in rapeseed/canola. Several disease screening methods including petiole inoculation technique (PIT) (Zhao et al. 2004;Bradley et al. 2006), detached leaf inoculation (Zhao and Meng 2003;Wu et al. 2013), detached or intact stem inoculation Wu et al. 2013;Wei et al. 2016;Qasim et al. 2020;Roy et al. 2021;Shahoveisi et al. 2021) have been used to evaluate the genetic resistance of rapeseed/ canola germplasm at different developmental stages under controlled and field environments. Moreover, rapeseed/canola cultivars differ in their plant architecture, growth habits, and maturity (Bradley et al. 2006;Arifuzzaman and Rahman 2020;Rahman et al. 2022). Cultivated canola express spring (no vernalization needed to induce flowering), semi-winter (shorter period of vernalization require to induce flowering), or winter (vernalization needed over the winter to induce flowering) growth habits (Wang et al. 2011;Arifuzzaman and Rahman 2020;Rahman et al. 2021). In North Dakota, the leading canola-producing state in the USA, only spring canola is cultivated due to the shorter growing season. Poor winter hardiness prevents the cultivation of semi-winter and winter ecotypes canola. It is difficult to use the stem inoculation method to screen all ecotypes of B. napus because of the challenge to synchronous inoculation time and the vernalization requirement. Therefore, a quick, efficient, and reliable early growth stage screening method is needed that would allow to circumvent bolting, flowering, or vernalization issues. This would ultimately accelerate the screening process and allow the simultaneous screening of genotypes with any of the three growth habit types. Here we implemented a PIT screening at four to five leaf stage seedlings to evaluate a diverse set of rapeseed/canola germplasm for resistance to S. sclerotiorum.
SSR resistance is a quantitatively inherited trait controlled by polygenes with minor additive and partially dominant effects, which are affected by the environment (Wei et al. 2016;Wu et al. 2016;Qasim et al. 2020;Roy et al. 2021;Derbyshire et al. 2021). Genetic mapping studies of SSR resistance through quantitative trait loci (QTL) analysis, based on classical linkage mapping strategy with bi-parental mapping populations, were commonly used for the purpose of identifying functional genes, and to position molecular DNA markers associated with SSR resistance. A majority of the SSR resistance QTL were located on chromosomes A01, A02, A03, A06, A7, A08, A09, C01, C02, C03, C04, C06, C08, and C09 (Zhao and Meng 2003;Zhao et al. 2006;Yin et al. 2010;Wu et al. 2013;Wei et al. 2014;Behla et al. 2017;Qasim et al. 2020;Shahoveisi et al. 2021) using a PIT, and detached leaf and/or stem inoculation technique at various developmental stages. Despite these efforts, no major QTL or gene conferring resistance to SSR was fine mapped or cloned, which seriously limits the research into genetic manipulation of SSR disease resistance breeding. Moreover, bi-parental mapping population lacks allelic diversity and has a fewer number of recombination events which limits mapping resolution (Korte and Farlow 2013). Genome-wide association (GWA) mapping has emerged as a robust approach to dissect complex traits and identify novel and superior alleles by capturing the resistance diversity in a germplasm collection to be utilized in marker-assisted breeding. GWA mapping is based on the utilization of linkage disequilibrium (LD) within a diverse population of genotypes which have undergone extensive historical and evolutionary recombination events leading to the development of shortened LD segments. Abundant genetic allelic diversity and faster LD decay provides more promising opportunities to achieve high mapping resolution for the significant SNP/ marker-trait associations (MTAs) than the traditional linkage mapping (Nordborg and Weigel 2008).
To date, only a few GWA studies have been carried out to identify MTAs for mapping SSR resistance in rapeseed/ canola (Gyawali et al. 2016;Wei et al. 2016;Wu et al. 2016;Roy et al. 2021). The screening procedure used in these studies were based on direct inoculation of mycelium to the main intact stem of the growing plant and/or to its detached stem during the flowering stage. Several investigations on SSR resistance reported a significant interaction and negative correlation between flowering time and stem resistance (Wu et al. 2019;Zhang et al. 2019;Roy et al. 2021). Therefore, we screened our association panel using a PIT method to evaluate the performance of the germplasm collections at the seedling stage to eliminate the conflicts between early maturation and SSR resistance in rapeseed/canola breeding. To the best of our knowledge, this report is the first GWA analysis to identify useful quantitative trait nucleotides (QTNs) associated with SSR resistance in rapeseed/canola using the PIT disease screening method at the seedling stage.
One main issue with association mapping (AM) is the low power of detecting rare variants with small effects, which might be associated with economically important traits (Bernardo 2016). Typically, individual SNPs/QTLs identified from GWA and bi-parental linkage mapping explain less than 12% of the phenotypic variance with some exceptions, which suggests uncaptured genetic potential for SSR resistance breeding in B. napus remains to be discovered. To obtain the maximum genetic potential for SSR resistance traits, genomic selection (GS) has emerged as a promising genomics-assisted breeding approach that accounts for both major and minor QTL effects into the prediction framework. GS utilizes the full genome information regardless of its significance, for genomic-enabled prediction of the superior genotypes as a candidate for selection (Meuwissen et al. 2001;Crossa et al. 2017). GS combines genome-wide molecular markers and phenotypic data of a training population to develop a statistical model that predicts the breeding and/or genetic values of selection individuals/candidates which are only genotyped (Meuwissen et al. 2001;Crossa et al. 2017;Derbyshire et al. 2021). Previous GS studies for various agronomic traits, including blackleg and S. sclerotiorum disease resistance, have shown the potential of genomic prediction to accelerate the rapeseed/canola breeding (Würschum et al. 2014;Fikere et al. 2020;Roy et al. 2021;Derbyshire et al. 2021). The predictive abilities on adult plant resistance against S. sclerotiorum in rapeseed/canola by Derbyshire et al. (2021) and Roy et al. (2021) clearly indicate the potential of GS for improving complex SSR resistance. This further motivated us to explore the effectiveness of GP in predicting SSR resistant genotypes using B. napus plants at the seedling stage.
In this work, we used a diverse panel of 337 B. napus lines with the aim to (i) identify new sources of SSR resistant genotypes at the seedling stage; (ii) detect significant genomic regions, SNPs/QTNs associated with SSR resistance at the seedling stage by performing single-locus and multi-locus GWA models; and (iii) to assess the potential of GP for seedling stage SSR resistance in rapeseed/ canola.

Germplasm collection
A total of 337 diverse B. napus germplasm accessions, and breeding lines with worldwide geographical origin of 23 countries, were obtained from the North Central Regional Plant Introduction Station (NCRPIS), Ames, Iowa, USA, and North Dakota State University (Suppl. Table S1). The plant materials consisted of spring, semi-winter, and winter ecotypes/growth habits of rapeseed/canola. The experiments were conducted in the Agricultural Experiment Station Research Greenhouse Complex, North Dakota State University, Fargo, ND, USA, during 2019 and 2020. Plants were grown in the greenhouse at 22 ± 2 ℃ temperature with a 16-h photoperiod provided by natural sunlight supplemented with 400 W HPS PL 2000 lights (P.L. Light Systems Inc.).

Experimental design, inoculum preparation and Sclerotinia sclerotiorum disease phenotyping
The experiments were conducted twice using a randomized complete block design (RCBD) with three replications in each experiment. The first and second experiments were carried out during the year of 2020 (E1) and 2021 (E2) greenhouse cycle. For each replicate, six individual plants were inoculated with S. sclerotiorum pathogen, which resulted in a total of 36 (6 plants × 3 replications × 2 experiments) plants being evaluated for each genotype. In each experiment, all genotypes were screened in three batches with two commercially available spring canola hybrid cultivars, 'Pioneer 45S51' and 'Pioneer 45S56' as resistant checks, and the publicly available Canadian cultivar 'Westar' as a susceptible check.
Throughout the study, a highly virulent single isolate WM031 of S. sclerotiorum was used for all inoculations to germplasm accessions Shahoveisi et al. 2021). The inoculum was prepared by culturing the surface-sterilized sclerotia of the isolate on autoclaved potato dextrose agar (PDA) medium (24 gL -1 potato dextrose broth and 15 gL-1 agar) at 22-24 ℃. Then mycelium plugs from the actively growing edges were sub-cultured on another PDA plates at room temperature for 48 h. Threeweek-old 4-5 leaf stage seedlings were inoculated using the PIT method described by Zhao et al. (2004) with a few modifications. In brief, the petioles of the second fully expanded leaf of the seedling were excised 2.5 cm from the main stem using scissors. Then two mycelium plugs of sclerotinia isolates from the actively growing edges of growing mycelium were loaded into a sterilized 200 µl pipette tip by pushing the open end of the micropipette tip into the 48-h-old culture plates. After that, the loaded tips were carefully pushed onto the severed petioles making sure the inner side of the agar plug flush with the top of the petiole tip (Fig. 1a). Two separate disease scoring systems were used to classify the phenotypic response. The inoculated plants were observed for two weeks, and the response of individual plants of each line was determined by days to wilt (DW). A plant was considered wilted when the infected main stem girdled completely or the leaves of the infected plant became irreversibly flaccid (Fig. 1f). DW were recorded daily for the following two weeks starting on the third day after inoculation. In addition to DW, lesion phenotypes (LP) were scored in a 1-to-5 rating scale at 3, 4, and 7 days post-inoculation (dpi) denoted as LP_3dpi, LP_4dpi, LP_7dpi. The phenotypic response was categorized according to Zhao et al. (2004) with following modifications: 1 = unaffected, no symptoms on the main stem; 2 = slightly affected, small size lesions (≤ 1.0 cm) at junction of petiole and stem, no water-soaked lesion, no wilt; 3 = moderately affected, small water-soaked lesions 1 3 (1.0-≤ 2.0 cm), no wilt; 4 = severely affected, expanded and sunken water-soaked lesion (≥ 2.0 cm), no wilt; and 5 = dead, expanded, sunken, and water-soaked lesion resulting complete wilting or topple over of the infected plant ( Fig. 1b-f).

Phenotypic data analyses
Data on DW and LP on 3 (LP_3dpi), 4 (LP_4dpi), and 7 (LP_7dpi) days were subjected to analysis of variance (ANOVA) for experiment-wise and combined across experiments in SAS version 9.4 (SAS Institute, Cary, NC). Data from the both experiments were combined if the ratio of the effective error variance for each trait was less than tenfold (Tabachnick and Fidell 2000;Rahman et al. 2019;Arifuzzaman and Rahman 2020). Best linear unbiased predictions (BLUPs) for all studied traits were used as the phenotypic values for the subsequent GWA and genomic prediction analyses. The BLUP was calculated using PROC MIXED procedure with the restricted maximum likelihood (REML) option in SAS 9.4. To obtain BLUP values combined across experiments, experiment was considered as fixed effect and genotype, genotype-by-experiment interaction, and replication nested within experiment was used as random. The broad-sense heritability (H 2 ) on family mean basis was calculated for each trait by estimating variance components using PROC MIXED procedure with REML option in SAS 9.4 (Holland et al. 2003). All effects were assumed to be random. Broad-sense heritability on an entry mean basis was calculated as follows: where σ 2 g , σ 2 ge , and σ 2 e represent the genotype, genotypeby-experiment interaction, and residual error variances, respectively; n and r were the number of experiments and replicates per experiments, respectively.
Descriptive statistics of raw values for SSR traits were computed using the PROC UNIVARIATAE procedure in SAS 9.4. Pearson's correlation was conducted to examine the relationship between all traits using R 4.1.2 (R Core Team 2021) using the raw phenotypic values. The correlation heatmap and density distribution plots of phenotypic value were performed using 'corrplot' and 'ggplot2' packages, respectively in R.

Genotyping data and quality control
Each member of the diversity panel was genotyped as described by Roy et al. (2021) and Rahman et al. (2022). In brief, total genomic DNA was extracted from fresh young and lyophilized leaf tissues using Qiagen DNeasy kit (Qiagen, CA, USA). The extracted DNA was quantified with a NanoDrop 2000/2000c Spectrophotometer (Thermo Fisher Scientific) and diluted to 50 ng/µl. Genomic libraries were prepared using ApekI enzyme digestion described by Elshire et al. (2011). The library was sequenced as single-end reads at the University of Texas Southwestern Medical Center, Dallas, Texas, USA, using Illumina Hi-Seq 2500 sequencer. Bowtie 2 (version 2.3.0) (Langmead and Salzberg 2012) was used to align the single-end sequencing reads against the 'ZS11' reference genome sequence (Sun et al. 2017).
n + 2 e nr TASSEL 5 GBSv2 pipeline (Glaubitz et al. 2014) was used to call the bi-allelic variant, which resulted in 497,336 unfiltered SNPs. Low-quality SNP markers were filtered with an individual read depth greater than 3, missing data less than 25%, minor allele frequency (MAF) greater than 5%, and physical distance (thin) less than 500 bp with VCFtools (version 0.1.16) (Danecek et al. 2011). Since canola is a self-pollinated crop, SNPs that were more than 25% heterozygous were removed using TASSEL 5.0 (Bradbury et al. 2007). After applying quality filtering, a total of 38,510 high-quality SNPs were obtained. The SNPs located outside of the chromosomes (unknown position) were removed, retaining 34,519 SNPs. SNP loci with missing values were imputed in Beagle 5.1 (Browning et al. 2018).

Single-locus genome-wide association analyses
SNPs with less than 5% MAF were removed from a total of 34,519 markers, leaving 27,282 high-quality SNPs for subsequent GWA analyses (available: https:// figsh are. com/s/ 018c2 34947 57b55 99b77). GEMMA software (version 0.98.1) (Zhou and Stephens 2012) was used for single-locus (SL) GWA analyses using a mixed linear model (MLM). The first three principal components (PCA) calculated by prcomp () function in R were embedded as covariates in the GWA analyses to control the confounding effect of population structure. Model-based clustering was performed using the first three PCAs to determine the subpopulations among the association panel with the Mclust package in R. The GEMMA-MLM was executed with the following command in the GEMMA (version 0.98.1) software: 'gemma -g [geno- .' A kinship matrix was incorporated as a random effect. The matrix was computed using the centered relatedness procedure in GEMMA. The significance threshold was determined using the method proposed by Li and Ji (2005) to determine the significant threshold value for the identified SNPs. In this method, for the 27,282 SNPs we calculated the effective number of independent loci (M eff ) by estimating correlation matrix and eigenvalue decomposition. The test criteria were then adjusted using the M eff with the following correction by Sidak (1967): where α p is the comparison-wise error rate and α e is the experiment-wise error rate (α e = 0.05).

Multi-locus genome-wide association analyses
Multi-locus (ML) GWA analyses were implemented using three multi-locus GWA algorithms that include MLMM p = 1 − e 1∕M eff 1 3 (Segura et al. 2012), FarmCPU , and mrMLM (Wang et al. 2016). For the ML models, we selected the three PCA to control population genetic stratification that we used for GEMMA-MLM. The MLMM and FarmCPU models were carried out using the GAPIT (version 3.0) R package (Wang and Zhang 2021). The mrMLM GWA model was implemented using the R package 'mrMLM' (Wang et al. 2016) with default parameters. The critical significant threshold between a trait and SNPs for all ML models was set to P ≤ 1.0 × 10 -3 [− log10 (P) ≥ 3.0], which has been broadly adopted by other researchers in various studies (Li et al. 2018;Xu et al. 2018;Karikari et al. 2020). The GWA results were visualized with Manhattan plot and comparative quantile-quantile (Q-Q) plots by plotting the observed P values against expected P values generated using the CMplot package in R language (https:// github. com/ YinLi Lin/R-CMplot).

Candidate gene search
The significant MTAs identified in at least two traits or two or more GWA models (either SL or ML) were selected for potential candidate gene search that may be associated with disease resistance. Candidate genes were searched within the LD blocks, where associated significant SNPs were located, and were regarded as the candidate gene search interval. If the detected SNPs were not located in the LD block, genomic regions spanning ± 50 kbp flanking regions of the significant MTAs were used as potential candidate gene interval in B. napus 'ZS11' reference gene models (Sun et al. 2017). LD block analyses on the same chromosome were computed by Haploview v4.1 with the default settings (Barrett et al. 2005).

Genomic prediction
Six genomic prediction models including rrBLUP, and five Bayesian models such as Bayes A (BA), Bayes B (BB) (Meuwissen et al. 2001), Bayes C (BC) ), Bayesian LASSO (BL) ), and Bayesian Ridge Regression (BRR) (Meuwissen et al. 2001), were used in this study. All models were performed in R language. The GP model rrBLUP was constructed using the package 'rrB-LUP' (Endelman 2011) and BGLR (version 4.0.4) package was used to fit the Bayesian GP models (Pérez and de los Campos 2014). All the analyses for the Bayesian models were performed for 5000 Monte Carlo Markov chain iterations with a 1,000 burn-iterations. The construction of genomic prediction models can be represented with the following general linear regression equation: where y is the vector of phenotypic values, µ is the intercept/ grand mean, X is the standardized marker genotype matrix, y = + X + β is the estimated random additive marker effects, and ε is the residual error term.
In brief, the rrBLUP model assumes that all the markers effects are normally distributed and all these marker effects have identical variance (Meuwissen et al. 2001), whereas Bayesian models may utilize different prior distributions which result in different levels of effect size shrinkage with various proportions of zero effect markers (Habier et al. 2011;Desta and Ortiz 2014;Crossa et al. 2017). In BA, the effect of each marker is estimated from a Gaussian distribution and markers are assumed to have different variances. The BB model is similar to BA, but allows some of the marker effects with zero variance. The BC model assumes a priori that markers have normally distributed effects with probability π and no effect with probability (1-π) (Meuwissen et al. 2017). The BL model applies both shrinkage and variable selection. The marker effects of the BL method are estimated from a double exponential distribution ). The BRR creates equal shrinkage of all the marker effects toward zero and produces a Gaussian distribution of the marker effects (Desta and Ortiz 2014).
The predictive ability of the GP models was tested with fivefold cross-validation (with 270 individuals as training set and remaining 67 individuals as validation set in each fold) and replicated 100 times to avoid biases in the estimation. Predictive ability of each trait is calculated as the Pearson correlation (r) between the average of the predicted genomic estimated breeding values (GEBVs) and the BLUPs in all the cross-validation sets. Prediction accuracy is defined as the correlation between GEBVs and the true breeding values. Since the true breeding values of tested traits are not known, we approximated the prediction accuracy by the correlation between the GEBVs and the observed phenotypic values (BLUP) divided by the square root of the phenotypic heritability ( √ H 2 ) following Lorenz et al. (2011) andJarquín et al. (2014). Thus, the accuracy of the models was estimated by dividing the mean of Pearson's r between GEBVs and phenotype values from all cross-validations with 100 cycles by square root of broad-sense heritability (H).

Phenotypic variations for Sclerotinia sclerotiorum reactions and correlation among phenotypic traits
SSR disease reactions can be variable under field environments; therefore, phenotyping the collection of germplasm against S. sclerotiorum was performed in the greenhouse under a controlled environment. A continuous and broad range of phenotypic variations were observed for days to wilt (DW) and lesion phenotypes (LP) traits among the genotypes in the study (Fig. 1a-f;  Fig. S1; Table 1). The BLUP values for DW varied from 4.0 to 8.9 days with an overall mean of 5.4 days and coefficient of variation is 25.3%. The variations of BLUP values observed for LP scores at 3, 4, and 7 dpi, ranged (mean) from 2.2 to 3.9 (2.8), 2.8 to 4.6 (3.8), and 3.9 to 4.9 (4.8), respectively (Table 1; Suppl. Table S2). The coefficient of variation (CV) of LP scores of the association population at different days varied from 6.4 to 27.2.% (Suppl . Table S2). Based on the phenotypic data, a few genotypes, which performed better than the resistant check cultivars used in this study, were identified as promising sources of resistance to SSR at the seedling stage. The BLUP values of the top five promising resistant genotypes ranged from 6.7 to 8.9 for DW, 2.2 to 2.4 for LP_3dpi, 2.9 to 3.0 for LP_4dpi, and 3.9 to 4.5 for LP_7dpi. However, the observed phenotypic responses of the resistant checks 'Pioneer 45S51' were 4.9, 3.2, 4.1, and 4.9 for DW, LP_3dpi, LP_4dpi, and LP_7dpi, respectively, and 'Pioneer 45S56' were 5.3, 2.9, 3.8, and 4.9 for DW, LP_3dpi, LP_4dpi, and LP_7dpi, respectively. The phenotypic response of susceptible check 'Westar' cultivar was 4.1, 3.8, 4.5, and 4.9 for DW, LP_3dpi, LP_4dpi, and LP_7dpi, respectively (Table 1; Suppl. Table S1). Analysis of variance (ANOVA) for SSR reaction in terms of DW and LP scores on different days revealed significant differences (P ≤ 0.001) among the genotypes, and interaction of genotype by experiment (Suppl . Table S3).
Highly significant correlations were observed among the phenotypic traits for SSR reaction. For instance, significant negative associations were found for DW with LP_3dpi (r = − 0.76), LP_4dpi (r = − 0.90), and LP_7dpi (r = − 0.83) at P ≤ 0.001 (Fig. S2). The estimated broadsense heritability of SSR resistance on entry mean basis across the two experiments was 0.69, 0.67, 0.68, 0.61 for DW, LP_3dpi, LP_4dpi, and LP_7dpi, respectively (Table 1). Medium to high heritability for SSR resistance in the phenotypic traits indicated that the phenotypic data were suitable for further genetic analyses.

SNP distribution and population structure analysis
After quality filtering and removal of markers with MAF < 5%, a total of 27,282 high-quality SNPs were used in the current study. These SNPs span a length of 854.3 Mb genome sequence representing 75.6% coverage of the B. napus genome (~ 1130 Mb). The number of SNPs were uneven among the 19 chromosome and ranged from 714-2386 SNPs per chromosome with the average SNP per chromosome was 1436, where the chromosomes 4 and 13 have the lowest (714 SNPs) and highest (2386 SNPs), respectively, while the average SNP per chromosome was 1436. The mean SNP density was approximately one SNP per 31.3 kb (Fig. 2a). Based on the 27,282 markers, principal component analysis (PCA) and kinship analyses were performed to identify the underlying genetic differences of the genotypes. The first three PCA explained 22.2% of the genotypic variation and were included in the GWA mapping model to control the confounding effect of population stratification. Furthermore, model-based clustering analysis using the first three PCA identified five subgroups within the genotypes based on three ecotypes (Fig. 2b).

Marker-trait associations detected for SSR resistance by multi-locus GWA analyses
Three multi-locus (ML) GWA algorithms: MLMM, Farm-CPU, and mrMLM, detected a total of 208 SNPs corresponding to 205 loci across all the 19 chromosomes of B. napus genome [− log10 (P) = 3.0 − 13.1] (Suppl. Table S4). The number of SNPs detected by the three ML-GWAS methods ranged from 10 to 46. The highest number of 46 SNPs was detected for LP_7dpi trait by MLMM, whereas the lowest number of 10 SNPs was found to be associated for LP_3dpi by mrMLM method. A total of 27 out of 208 SNPs were identified simultaneously in at least two phenotyped SSR resistance traits by two or more ML methods for any one of the traits. The estimated allelic effects ranged between − 0.36 to 0.40, − 0.18 to 0.14, − 0.15 to 0.13, and − 0.12 to 0.08 for DW, LP_3dpi, LP_4dpi, and LP_7dpi traits, respectively. The explained phenotypic variation accounted for by the significant SNPs ranged from 2.2-6.0%, 2.72-9.72%, 1.25-12.22%, and 2.65-10.13% for DW, LP_3dpi, LP_4dpi, and LP_7dpi traits, respectively (Suppl.

Commonly identified marker-trait associations among the SSR resistance traits, among and between single-locus and multi-locus GWAS methods
Of the 41 detected SNPs by SL GEMMA-MLM method, there were a maximum of 21 SNPs for LP_7dpi, 12 SNPs for LP_4dpi, 10 SNPs for DW, and 9 SNPs for LP_3dpi traits (Fig. S3). Five SNPs were mutually identified between DW and LP_4dpi, followed by 4 SNPs between DW and LP_7dpi, 3 SNPs between LP_3dpi and LP_4dpi trait, and only single SNP between DW and LP_3dpi (Suppl .  Table S4, S5). Moreover, one SNP SCM002771.2_77997199 on chromosome C03 for DW, LP_3dpi, and LP_4dpi traits and another SNP SCM002776.2_3494019 on chromosome C08 for DW, LP_4dpi, and LP_7dpi traits were found to be co-localized by the SL method. All of the QTNs detected with the SL method were also associated with the four SSR resistance traits by any of the ML-GWA models. In addition to the 41 QTNs identified by SL, GWA analyses by ML methods detected additional 167 SNPs associated with SSR phenotypic traits. The number of identified QTNs by all the ML models for SSR resistance traits varied between 63 and 70, whereas the number of QTNs for each of the ML models ranged between 10 and 46. The highest (46) number of QTNs was detected by MLMM(LP_7dpi), and the lowest (10) QTNs by mrMLM model out of the three ML models for LP_7dpi (Fig. S3). Comparison of the ML models demonstrated that each model has the power to detect QTNs concurrently from each other and a few QTNs (ranged 2 to 9) were detected by all the three models for each trait. However, to obtain more reliable results, only the SNPs that were simultaneously detected by a combination of SL and any of the ML methods or at least two of the ML methods or at least two traits were considered as significant QTNs. Thus, a total of 98 QTNs controlling SSR resistance traits were obtained (Suppl. Table S5). These QTNs will serve as a valuable source and could provide promising opportunities to facilitate cost-effective MAS breeding of SSR resistance. Manhattan and Q-Q plots summarizing the GWA results of all the phenotypic traits for SSR resistance by SL (GEMMA-MLM) and ML (MLMM, FarmCPU, mrMLM) algorithms are present in Fig. 3a-d and Figs. S4, S5, S6, and S7. All GWA models were compared with the studied phenotypic traits to determine whether the models control false positives and false negatives. The Q-Q plot depicts the expected negative log10 (P) values versus the expected negative log10 (P) values across all markers. Q-Q plots of models including GEMMA-MLM and MLMM had a straight line with slightly deviated tail, which indicated that these two models reduced false positives (Fig. 3a-c). However, most of the SNPs were close to the straight line or a little bit inflated downward, indicating that they might have been reported as false negatives (Fig. 3a-d). In contrast, examination of Q-Q plots of FarmCPU and mrMLM models showed a sharp upward deviation from the expected P value distribution in the tail area, indicating these models controlled both false positives and false negatives (Fig. 3a-c).

Candidate gene prediction
To identify the potential candidate genes for the SSR resistance, the significant SNPs detected in at least two traits or with two or more GWA models were used for candidate gene mining using the 'ZS11' reference genome sequence database (Sun et al. 2017). With this criterion, 83 putative candidate genes with known functions associated with plant disease resistance mechanisms were identified. To discover putative biological functions, the annotated protein of the corresponding candidate gene model were searched on Uniport database (https:// www. unipr ot. org/ unipr ot/) (Suppl. Table S6). The biological processes of the detected candidate genes were involved in defense response, defense response to fungus, response to a molecule of fungal origin, response to chitin, programmed cell death, callose deposition in cell wall, response to salicylic acid, indole glucosinolate biosynthetic process, induced systemic resistance, jasmonic 1 3

Fig. 3 Circular Manhattan plots
showing statistically significant SNPs based on single-locus (SL) GEMMA-MLM, and three multi-locus (ML) MLMM, FarmCPU, and mrMLM models located on 19 chromosomes for Sclerotinia sclerotiorum resistance at the seedling stage. Associations for the days to wilt (a), lesion phenotype (LP) scores at 3 days post-inoculation (dpi) (b), LP_4dpi (c), and LP_7dpi (d) are shown. A multi-track Q-Q plot for each trait with all the GWA models is presented to the right of each Manhattan plot. The threshold values for SL and ML models were set up at − log 10 (P) ≥ 3.6 and − log 10 (P) ≥ 3.0, respectively acid-mediated signaling pathway, ethylene-dependent systemic resistance, systemic acquired resistance, pattern recognition receptor signaling pathway, response to wounding, protein kinase activity, response to oxidative stress, toxin catabolic process, immune response, reactive oxygen species metabolic process, brassinosteroid-mediated signaling pathway, and other biological processes which might play a key role in early-stage SSR resistance in rapeseed/canola (Suppl . Table S6).

Genomic prediction (GP)
The mean predictive ability and prediction accuracy of six GS models are shown in Fig. 4a-

Discussion
In the present study, the reaction of 337 rapeseed/canola genotypes against S. sclerotiorum was evaluated for DW and LP at 3, 4, and 7dpi to identify the potential SSR resistant genotypes. The LP scores were recorded at different time points to distinguish differences in disease progress among the genotypes and to identify a single time point for phenotypic scoring that correlates with the commonly used DW data for SSR phenotyping. QTL mapping studies conducted by Zhao et al. (2006) used DW and stem lesion length data at 4 dpi, whereas Behla et al. (2017) used only DW. The use of multiple phenotyping scoring systems would provide valuable insights to accurately evaluate the resistance performance of the genotypes and detect additional QTNs associated with disease resistance Shahoveisi et al. 2021). Furthermore, the evaluation of extensive phenotyping with multiple scoring systems would enable researchers to select a single time point to score disease phenotype. Based on our phenotyping screening results, a wide range of phenotypic variability were recorded in response to the S. sclerotiorum infection in the studied rapeseed/canola germplasm. The continuous and broad range of observed phenotypic responses reinforce the notion that SSR resistance in B. napus is a quantitatively inherited trait, controlled by many minor genes with small effect (Zhao et al. 2006;Wu et al. 2013Wu et al. , 2016Wei et al. 2016;Qasim et al. 2020;Roy et al. 2021). Pearson correlation analyses among the phenotypic traits revealed that DW trait was significantly and strongly correlated with the different time points LP scores data (r = − 0.76 to − 0.90). However, LP score at 4 dpi was found to have the strongest association (r = − 0.90) with the DW. The strong association among DW and LP scores was also reported by Zhao et al. (2004Zhao et al. ( , 2006. Therefore, based on these findings, LP_4dpi could be used as a proxy criterion to DW when evaluating rapeseed/canola germplasm for resistance to S. sclerotiorum. The estimated broad-sense heritability (H 2 ) for SSR resistance in the 337 B. napus germplasm (61-69%) implied the majority of the observed phenotypic variation was controlled by genetic factors. This level of heritability is consistent with the previous SSR resistance studies (Zhao et al. 2006;Wu et al. 2013;Wei et al. 2016;Qasim et al. 2020;Roy et al. 2021), which further indicate that phenotypic selection is effective and therefore suitable for subsequent GWA analyses to detect favorable alleles conferring SSR resistance to facilitate MAS breeding. The power of detecting MTAs by GWA study is limited by several factors including population size, cryptic population structure, linkage disequilibrium, heritability of the trait, underlying genetic architecture of the trait of interest, and the statistical models used (Gupta et al. 2005;Yu et al. 2006;Josephs et al. 2017). The diversity panel used here consists of 337 rapeseed/canola genotypes originating from 23 countries and comprising 0.3% Australian, 28.8% Asian, 35.0% European, and 35.9% North American origin. In terms of growth habits, germplasm accessions were comprised of 42.7% spring, 16.9% semi-winter, 39.5% winter, and 0.9% of unknown ecotypes. Therefore, the geographical distribution of the used genotypes in our study provides an ideal diverse panel with all forms of ecotypes collected from the major rapeseed/canola growing regions. The high phenotypic variability among the genotypes coupled with an ideal diversity panel with good worldwide geographical coverage, high mean SNP density (~ one SNP per 31.3 kb) enhances QTN detection via GWA analyses. Hereafter, we implemented the first GWA study at the seedling stage to identify the significant SNPs, genomic regions, and putative 1 3 Fig. 4 Violin plot with predictive ability and accuracy of six genomic prediction models, rrBLUP, Bayes A (BA), Bayes B (BB), Bayes C (BC), Bayesian Lasso (BL), and Bayesisan ridge regression (BRR) to detect sclerotinia stem rot resistant genotypes using phenotypic data of days to wilt (a), lesion phenotypes (LP) at 3 days postinoculation (dpi) (b), LP_4dpi (c), and LP_7dpi (d) obtained from canola/rapeseed plants inoculated at the seedling stage with S. sclerotiorum using the petiole inoculation method. The red triangle in each violin plot represents the mean predictive ability. The numbers above red triangle represent the predictive ability (r) at the top and prediction accuracy in brackets (colour figure online) candidate genes conferring SSR resistance in rapeseed/ canola.
As the pattern of quantitatively inherited SSR resistance is complex and controlled by many genes with small effect, we simultaneously used multiple GWA models, i.e., one SL (GEMMA-MLM) and three ML (MLMM, FarmCPU, mrMLM) GWA models to discover SNPs associated with the traits. Although SL model (MLM) is widely used to detect the genetic variants for traits of interest in many crop species, it has several limitations to dissect complex traits. SL models perform one-dimensional genomic scan by testing one marker at a time and also fail to simultaneously match the true overall genetic model of complex/quantitative traits controlled by multiple loci. False negative (Type II error) could also result from the MLM-based SL models due to the model overfitting, where some potentially important associations could be missed ). On the other hand, tremendous statistical improvement efforts were made over the few years to overcome the problems associated with SL GWA models for the dissection of complex traits. Several multi-locus (ML) GWA algorithms, such as MLMM (Segura et al. 2012), FarmCPU , mrMLM (Wang et al. 2016), FASTmrEMMA (Wen et al. 2018), ISIS EM-BLASSO (Tamba et al. 2017), and pLARmEB , were developed. The advantages of these ML methods are that no multiple test correction is required due to the multi-locus nature of the model, and also have more statistical power and accuracy to detect associations than SL models (Wang et al. 2016;Xu et al. 2017). Similar trends were also observed in the current study, where for most of the traits, ML-GWA models detected more QTNs compared to the SL model. Among all the GWA models implemented, MLMM and FarmCPU have shown the highest power to detect MTAs than mrMLM and GEMMA-MLM.
GWA analyses revealed a total of 208 significant QTNs corresponding to 205 loci for all the studied SSR resistance traits. However, to obtain more reliable results, only SNPs simultaneously detected by both SL and any of the ML methods or at least two of the ML methods or two phenotypic traits were considered as significant QTNs. Thus, a total of 98 QTNs controlling SSR resistance traits were obtained. Additionally, 44 QTNs were simultaneously mapped for at least two SSR resistance traits. Use of SL and multiple ML models for GWA analyses with four SSR resistance phenotypic traits improved the reliability of QTNs detection and was also complementary to each other in identifying common and more significant QTNs for the trait of interest. The ML models detected more significant QTNs over SL models which confirm the power and robustness of ML-GWAS models compared to SL. Similar trends were also observed in cotton (Li et al. 2018), maize (Xu et al. 2018), and soybean (Kaler et al. 2020;Karikari et al. 2020), where ML models detected more significant SNPs over SL model. Therefore, implementing GWA analyses using ML models in conjunction with SL model would enable the detection of more QTNs associated with traits of interest and provide promising opportunities to facilitate the genomics-assisted SSR resistance rapeseed/canola breeding.
Presently, several QTLs/markers associated with early (seedling) and adult stage SSR resistance either in the form of stem, leaf, and days to wilt resistance have been identified using bi-parental linkage mapping and association mapping studies (Zhao and Meng 2003;Zhao et al. 2006;Yin et al. 2010;Wu et al. 2013Wu et al. , 2016Wei et al. 2014Wei et al. , 2016Gyawali et al. 2016;Zhang et al. 2019;Qasim et al. 2020;Roy et al. 2021;Shahoveisi et al. 2021). The physical position of the previously reported QTLs/MTAs was based on the 'Darmorbzh' reference genome (Chalhoub et al. 2014). However, in our current study, we used 'ZS11' as reference genome sequence (Sun et al. 2017), which was aligned with 'Darmorbzh' reference genome. The physical location comparison of our identified 208 SNPs with the previously reported QTLs/ markers revealed a total of 54 SNPs corresponding to previously reported SNPs/QTLs detected based on linkage and/ or association studies. Markers that are linked with QTLs for SSR resistance located within ~ 500 kb of the same genomic regions were considered as the same loci. In our current study, a total of 15 SNPs were detected on chromosome A09 with multiple GWA models and traits. Among them 2 SNPs (SCM002767.2_21520686, and SCM002767.2_21528172) (21.52-21.53 Mb), located ~ 7.5 kb apart, were in the close proximity of the mapped QTL by Wu et al. (2013), andQasim et al. (2020) for stem resistance at adult stage. Another SNP (SCM002767.2_27713481) on chromosome A09 (27.71 Mb) co-localized with the qSR11-1 QTL between the physical position of 27.13-29.36 Mb (Wei et al. 2014). Li et al. (2015) defined the genomic regions spanned 22.5-27.5 Mb on chromosome A09 as the conserved QTL regions for S. sclerotiorum resistance based on the integrated and comparative QTL analyses of the previously identified QTLs. Our findings also corroborate that this physical interval is rich in QTLs conferring resistance to SSR at both the seedling and adult plant stages. Moreover, these genomic regions might be a potential hotspot of stable QTL to carry out fine mapping or map-based cloning to take the full advantage of MAS in SSR resistance breeding. Two additional SNPs identified by multiple models and traits located on C02 (0.96-6.16 Mb) were also detected by other researchers. For instance, Zhao et al. (2006) detected Sll 12 (stem lesion length) and Dw 12 (days to wilt) QTL within the physical interval of 0.31-6.71 Mb for SSR resistance using PIT. More QTLs, such as qSR10-3 (1.03-3.95 Mb) and qSR11-2 (1.03-3.95 Mb) by Wei et al. (2014); qSRC2 (0.02-4.33 Mb) by Wu et al. (2019); and SRC2a (0.23-5.55 Mb) and SRC2b (0.12-5.19 Mb) by Qasim et al. (2020), were previously detected using stem inoculation 1 3 method as a screening technique for adult plant resistance. A GWA study conducted by Roy et al. (2021) identified 10 significant SNPs on A09 at 35.6-45.8 Mb for SSR resistance using the stem inoculation technique under field environments. Our GWA results also mapped 5 significant SNPs in this genomic region. The simultaneous detection of common QTNs controlling seedling and adult plant stage SSR resistance in different populations and genetic mapping methods provides valuable insights that warrants further exploration of these genomic regions to develop functional molecular markers that could potentially be used in MAS of target traits at all development stages. Besides the identification of the previously detected SNPs/QTLs, our study also revealed new genomic regions [A01 (1.58-12.96 Mb); A06 (2.50-7.47 Mb); C01 (11.45-16.51 Mb); C02  Mb)] that may contribute to better understanding of the architecture of S. sclerotiorum reaction and could provide more opportunities for SSR resistance breeding in rapeseed/canola.
The identification of stable QTNs/genomic regions is necessary to provide useful information to facilitate MAS. Therefore, the QTNs detected in two or more traits and two or more multiple GWA mapping algorithms were selected for mining potential disease resistance candidate genes for further gene cloning and functional verifications. A total of 83 defense response candidate genes surrounding the stable QTNs were identified using the 'ZS11' reference genome. These candidate genes were categorized on the basis of their functional characteristics for disease resistance from several databases. A TIR-NB-LRR receptor-like protein (Gene ID: LOC106378095) annotated as disease resistance protein RPP1 was located on chromosome C02 (25.89-26.03 Mb) for multiple traits with multiple GWA algorithms. This candidate was also detected as a potential candidate gene for adult plant stage SSR resistance under field environments by Roy et al. (2021). Recent findings on the lifestyle of S. sclerotiorum revealed that there is a brief biotrophic phase followed by a necrotrophic phase, which suggests S. sclerotiorum is a hemi-biotrophic pathogen rather than necrotrophic (Kabbage et al. 2015;Chittem et al. 2020). TIR-NB-LRR proteins regulate the activation of salicylic acid (SA)-dependent pathway to confer defense response against biotrophic and hemi-biotrophic pathogens. Therefore, SA may plays a positive role in the defense responses against S. sclerotiorum in rapeseed/canola, which agree with Weaver et al. (2006), and Nováková et al. (2014) findings. Another significant SNPs (SCM002761.2_479878) on chromosome A02 identified by multiple traits and multiple models are found to be linked with a germin-like protein subfamily 3 member 3 (GPL3) gene (Gene ID: LOC106434866). GLPs have diverse functions including contribution to plant defense reactions against different pathogens. Expression of GLPs was found to be associated with increased resistance against Sclerotinia sp. pathogen (Dong et al. 2008). Enhanced resistance in Arabidopsis thaliana plants against Verticillium longisporum and Rhizoctonia solani were also reported through the transgenic expression of BvGLP-1 by Knecht et al. (2010), whereas increased susceptibility was reported to the rice blast fungal pathogen when GLP genes were silenced (Manosalva et al. 2009). However, in B. napus BnGLP3 and BnGLP12 were found to be upregulated and generated more H 2 O 2 formation in the partially resistant cultivars compared to the susceptible cultivars upon S. sclerotiorum infection. The tolerant lines also generate an increase in H 2 O 2 leading to the oxidative burst at the early state of S. sclerotiorum-infected leaves and thereby restricting lesion formation compared with the susceptible cultivar (Rietz et al. 2012). Other candidate genes for SSR resistance would also provide useful insights for efforts to achieve S. sclerotiorum resistance rapeseed/canola germplasm.
The advantage of GS over GWA study is that it simultaneously exploits the predictive power of all the genome-wide distributed markers (Meuwissen et al. 2001). Therefore, GS is considered as an effective genomic strategy for the improvement of complex traits in crops that could potentially capture minor-to medium effect loci (Meuwissen et al. 2001;Würschum et al. 2014;Crossa et al. 2017). Since SSR resistance is a polygenic trait, controlled by numerous minor-effect QTL, MAS may not always be efficient. Therefore, GS could provide an alternative breeding technique more capable of accelerating breeding progress by increasing genetic gains. This led us to assess the predictive ability and accuracy for seedling stage SSR resistance by implanting six GS models. The predictive ability was moderate to high, depending on trait specifications. These results clearly demonstrated that genome-wide markers are efficient in predicting resistance of rapeseed/canola genotypes in response to S. sclerotiorum attack. The highest mean predictive ability 0.66-0.67 was observed for LP_3pi traits followed by LP_4dpi (0.62-0.64), DW (0.60-0.61), and LP_7dpi traits (0.42-0.44). The resulting predictive ability from our study was comparable or higher than that of Wu et al. (2016), Derbyshire et al. (2021), and Roy et al. (2021) reporting predictive ability for adult plant resistance. The differences in predictive ability could be attributed due to the differences in size and genetic diversity of the populations, molecular markers, extent of linkage disequilibrium, and the heritability of the traits (Daetwyler et al. 2010;Isidro et al. 2015;Crossa et al. 2017). Some traits (LP_3dpi, LP_4dpi, and DW) were predicted with relatively high predictive ability (0.61-0.67). Similar trends were observed by Roy et al. (2021), who reported 28-56% increase in predictive ability depending on the assessment of various SSR phenotypic traits. These results suggest that the predictive ability of SSR resistance could be dependent on the use of the phenotypic scoring systems at both seedling and reproductive stages. There was a little difference in predictive ability existed among the different GS models across traits, and no one models outperformed others. Similar observations were also reported for S. sclerotiorum resistance in rapeseed/canola Derbyshire et al. 2021) and soybean (de Azevedo Peixoto et al. 2017). The slightly observed differences in predictive ability are likely due to the models underlying assumptions. The result of this study clearly indicates that GS holds promising potential for improving SSR resistance by predicting genotypes as a potential donor/ parent at the early stage (seedling) in efforts to breed SSR resistant rapeseed/canola cultivars. Knowledge of predictive ability obtained from the evaluation of diverse B. napus germplasm collections would greatly contribute to breeding pipelines and an important initial step for the establishment of genomic selection (Jarquin et al. 2016;Crossa et al. 2017;Rolling et al. 2020). However, evaluation of the genotypes and testing of the GP models in the field screenings would definitely help to confirm the application of GS for SSR resistance as the performance of the genotypes would be more representative in the farmer's field.
Finally, the use of PIT as an inoculation method to screen diverse B. napus accessions consisting of three ecotypes, i.e., spring, semi-winter, and winter, against S. sclerotiorum attack enabled us to identify promising resistant genotypes at the early developmental stage. The identified resistant lines will provide a valuable source for canola breeding efforts to improve durable resistance by developing SSR resistant cultivars for the growers. Moreover, screening of genotypes at the same growth stage will help to eliminate/overcome the effects of several physiological traits, which were known to have an indirect effect for SSR resistance evaluation. In addition to common use of DW as a phenotypic response, the highly correlated single time point phenotypic score LP_4dpi could be used as a potential phenotypic trait for large scale phenotyping of SSR resistant genotypes. Due to the quantitatively inherited nature of SSR resistance, it is more likely that QTNs identified here had small effects on resistance in B. napus. The detected significant markers identified by GWA or QTL mapping could be converted into kompetitive allele-specific PCR (KASP) markers for SNP validation (Semagn et al. 2014). Our future efforts would be directed toward validating the effects of these QTNs and the development of tightly linked markers to facilitate the costeffective MAS resistance breeding in rapeseed/canola. Considering the potential of GS for the improvement of polygenic traits, this is the first report of exploring the feasibility of GP for seedling stage SSR resistance. High predictive ability coupled with high prediction accuracy demonstrates the promise of GS for improving resistance toward S. sclerotiorum in B. napus and will assist in the establishment of GS research by enhancing germplasm for disease resistance breeding. However, assessment of disease resistance and GS models under multilocation field environment might be warranted in the future for effective SSR resistance genomicsassisted rapeseed/canola breeding.