Two extreme parent lines with resistance J964 and susceptibility J902 were used to develop the double haploid (DH) population for this study. The cross was made between two parents, and F₁ was developed. A DH population was developed using anther culture of the F₁ population. A total of 181 DH lines were developed using anther culture and were used for trait analysis and QTL mapping. All materials were grown and phenotypic measurement were recorded on the experimental farm of Huazhong Agriculture University, Wuhan, China, in the winter-type oilseed rape growing season in 2015–2016, 2016–2017, 2017–2018 and for measuring flowering time of DH lines used for RNA-seq grown again in 2018–2019. Parents along with the DH lines were grown in three replications following randomized complete block design (RCBD). The distance between the plants was kept 21 cm in each row and 30 cm between the rows. The field was managed following normal breeding practices.

Parents along with DH lines were artificially inoculated with S. sclerotiorum for three consecutive years. A mycelial agar plug (7 mm in diameter) punched from the margin of a 2-day-old culture of S. sclerotiorum grown on potato dextrose agar (PDA) was used as the inoculum. Lesion length was noted in centimeters after 7 days of inoculation.

Flowering data were noted in 2017 for three replications of the DH population. When the 25% of plants in each replication of each genotype started flowering, it was noted as initial flowering. There was only a small gap in days from initial (25%) to final (100%) flowering, as most of the plants from the same lines start flowering at the same time. The initial flowering time was then used for the correlation between flowering time and lesion length. Flowering time was also noted in 2019 for selected plants that were later used for RNA-seq analysis. Stem width data were taken in 2018 from 30 random genotypes, including three replications for each genotype. The stem width (diameter) was measured using a Vernier caliper in millimeters from the same internode and slightly above the point where inoculant was applied. Regression analysis (simple linear regression model) was performed using Microsoft Excel with lesion length (as the dependent variable) to flowering time (as the independent variable) and was the same for plant width, with width replacing flowering time.

Genomic DNA was extracted from 181 DH lines and parents grown in 2016 using the CTAB method as described earlier (Murray and Thompson, 1980). A total of 181 genotypes selected were used for genotyping along with parents. Genotyping was performed using Brassica 6 K Illumina Infinium HD Assay SNP arrays (Illumina Inc., San Diego, CA, United States) according to the manufacturer’s protocol (Infinium HD Assay Ultra Protocol Guide¹) (Mason et al., 2017). Illumina GenomeStudio software (Illumina Inc., San Diego, CA, United States) was used to cluster and visualize all SNP array data (Supplementary Table S1) for further analysis. Each SNP was rechecked manually to see if any error was there during the clustering analysis.

The SNP array data processing for the DH population followed a bi-filtering method described by Cai et al. (2015) was used, and the markers with same genotypes for the DH population were considered redundant SNPs and were removed. We only used the non-redundant SNPs (nrSNPs) for linkage map construction with MAPMAKER/EXP 3.0 and MSTmap software as described by Cai et al. (2014). The genetic distances between the markers was calculated using kosambi mapping function. QTL analysis was conducted using composite interval mapping (CIM) (Zeng, 1994) with Windows version of QTL cartographer 2.5 software² with the same protocol described earlier (Wu et al., 2013). A significance threshold (P = 0.05) for QTL was determined through permutation analysis using 1000 repetitions. Forward and reverse regression analysis was applied for QTL detection. The confidence interval of QTL was measured at 1-LOD intervals surrounding the QTL peak.

Comparative RNA-seq was performed for the two lines, one resistant and the other susceptible, with a consistent pattern of lesions in all 3 years of artificial inoculation. The main stems of the plants were inoculated at three sites on the neighboring internodes 40–60 cm above the ground using mycelial agar plugs. Mock-inoculated plants were inoculated with agar plugs only without any fungal growth. Epidermal stem tissues were harvested from 10 mm beyond the inoculation site and 1 mm deep at 24, 48, and 96 hours post-inoculation (hpi). Samples from each plant were mixed together from all three sites and considered as one replication, and three plants were used for three replications at one time point from each line for disease inoculation. One plant was mock-inoculated in each line at each time point, making three replicates covering three time points for mock-inoculated resistant and susceptible lines. The sampled tissues were immediately frozen in liquid nitrogen and later stored at −80°C before use for RNA extraction.

RNA was extracted from all 24 samples by using an NEB extraction kit following protocols by the manufacturers. RNA degradation and contamination were checked on 1% agarose gels. RNA purity was tested using a NanoPhotometer^® spectrophotometer (IMPLEN, CA, United States). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States). Sequencing libraries were created using the NEBNext^® UltraTM RNA Library Prep Kit for Illumina^® (NEB, United States) following the manufacturer’s recommendations, and index codes were added to attribute sequences to each sample. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Novaseq platform, and 150-bp paired-end reads were produced.

Raw data (raw reads) in fastq format were first processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N and low-quality reads from raw data. At the same time, Q20, Q30 and GC content of the clean data were calculated. All downstream analyses were conducted with clean data with high quality. ‘Darmor-bzh’ reference genome and gene model annotation files (Chalhoub et al., 2014) were downloaded from the genome website. The index of the reference genome was built using Hisat2 v2.0.5, and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5³. We selected Hisat2 as the mapping tool because Hisat2 can generate a database of splice junctions based on the gene model annotation file and thus a better mapping result than other non-splice mapping tools.

Feature Counts v1.5.0-p3 (Liao et al., 2014) was used to count the read numbers mapped to each gene, and then fragments per kilobase of transcript sequence per million (FPKM) fragments mapped of each gene was calculated based on the length of the gene and read count mapped to this gene. The FPKM considers the effect of sequencing depth and gene length for the read count at the same time and is currently the most commonly used method for estimating gene expression levels. Differential expression analysis of two groups (including biological replicates per group) was performed using the DESeq2 R package (1.16.1) (Love et al., 2014). DESeq2 provides statistical routines for measuring differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using Benjamini and Hochberg’s method (Haynes, 2013) for monitoring the false discovery rate. Genes with an adjusted P-value < 0.05 found by DESeq2 were considered as DEGs.

The expression pattern of the genes obtained from RNA-seq analysis was validated by performing the Quantitative real-time polymerase chain reaction (qRT-PCR). The expression pattern of the nine selected genes was studied with qRT-PCR analysis using the same samples used for RNA-seq analysis. The qRT-PCR assay was carried out as explained earlier by Shahid et al. (2019) using three biological replicates repeated twice as technical replicates. SYBR-Green MasterMix was used to perform qRT-PCR experiment. The gene copy specific primers were designed using Premier-5 program. The relative cycle threshold (Ct) value was used to measure relative expression by the 2^{–Δ Ct} method. BnaA09g14410D (PP2A-1) gene was used as control to normalize the expression data (Zhao et al., 2019). The detail of genes and primers are listed in Supplementary Table S4.

Parent lines showed clear differences in lesion length, noted 7 days post-inoculation (dpi) in all 3 years (Figure 1). The DH population developed from F₁ of two extreme parents assayed after 7 dpi revealed continuous segregation. It can be observed from the lesion length distribution that this trait is controlled quantitatively. Disease lesion length distribution among the population recorded at 7 dpi for the 3 years ranged from 3.43 to 17.80 cm in extreme lines. The frequency distribution of lesion length in 3 years among the DH population is shown in Figure 2.

The correlation among the 3 years of data was noted to see the consistency of the resistance pattern in DH lines. The correlation was highly significant, indicating the same pattern of lesion development among the genotypes and confirming the credibility of the method used (Table 1). The heritability for SSR resistance is as high as 68.36% (Table 2), indicating that the phenotype is primarily attributable to genetic factors with less environmental effects.

There are evidences that flowering time in B. napus affects SSR resistance (Wei et al., 2014; Wu et al., 2019); therefore, we noted differences in flowering initiation and flowering time in our study to determine whether flowering time affects disease resistance. Flowering time was found to be significantly negatively associated with disease lesion length (Figure 3B). The coefficient of correlation was as high as −0.53, suggesting that the plants with early flowering have greater susceptibility to disease with increased lesion length.

Using the 1401 SNP polymorphic markers for the parents and DH population of 181 lines, a linkage map was constructed that covered a distance of 972 cM with the kosambi function (Table 3). Seventeen QTLs were identified in three seasons using a linkage map and phenotype data of 3 years for disease lesion length. All 17 identified QTLs are present on six chromosomes out of 19 chromosomes in B. napus.

Seven QTLs were identified in 2016 on chromosomes A2, A9, and C3, five QTLs in 2017 on chromosomes A2, C2, C3, and C4, and five QTLs in 2018 on chromosomes A9, C2, and C6 (Table 4). None of the QTLs was identified in all 3 years. The chromosomes having QTLs for 2 years were A9 (2016, 2018), A2 (2016, 2017), C2 (2017, 2018), and C3 (2016, 2017). Among these QTLs identified, SRA9a in 2016 lies between the confidence interval for SRA9a in 2018 with R² 14 and 11%, respectively. SRC3a for 2016 lies within the confidence interval of SRC3a in 2017 with R² 7 and 10%, respectively. SRC2a in 2017 also lies within the confidence interval of QTL SRC2a for 2018, with R² 7.73 and 6.81%, respectively. All of the pairs of common QTLs with overlapping confidence intervals showed additive effects from the same common parent. These overlapping QTLs in 2 years are assigned the same name. It can be concluded that the same genes may be involved in resistance in different growing seasons, which warrants further research.

For more confirmation about the QTL regions possessing the genes contributing to SSR resistance, we collected some evidence from some past studies and compared them with our findings (Figure 4). We collected information on QTLs and their physical positions in the linkage group based on marker information from past studies and compared it with our findings (Supplementary Table S2). The overlapping regions among the QTLs identified in this study and some previous studies were observed, and 4 out of six chromosomes having QTLs from this study showed overlapping regions with past studies. QTL SRC2a on chromosome C2 identified in two seasons (2017 and 2018) in the present study have shown overlapping confidence intervals with the QTLs identified by Zhao et al. (2006); Wei et al. (2014), and Wu et al. (2019). QTL SRC3a sharing overlapping confidence intervals in this study for two seasons also shows overlapping confidence intervals with the previous finding of Wu et al. (2019). The QTL SRA9a, which is also present in two seasons, has shown overlapping confidence intervals with QTLs identified by Wu et al. (2013). The QTLs identified only in one season on A9 SRA9b have an overlapping confidence interval with the QTL identified by Wei et al. (2014), and SRA2c has an overlapping confidence interval with Wu et al. (2013). The QTL SRC6 identified on chromosome C6 in 2018 showed an overlapping confidence interval with Zhao et al. (2006) and Wu et al. (2013). To identify the genes responsible for resistance to SSR, we conducted RNA-seq analysis.

Also, the QTLs identified in other studies were on chromosomes A01, A03, A06, A08, A10, C03, C08, and C09 (Zhao et al., 2006; Wu et al., 2013, 2016a, 2019; Wei et al., 2014, 2016).

On the basis of 3 years of lesion length data, two consistent lines of resistant and susceptible genotypes were selected for RNA-seq analysis. Both lines had clear differences in their lesion length patterns in all 3-year extremes ranging from 4.7 to 15.2 (Figure 5B). The lines selected also had similar flowering times, as we noted that flowering time also affected plant disease resistance patterns (Figure 5A). In this study, we only attempted to focus on the mechanisms directly affecting SSR resistance in B. napus, as both lines were at the same growth period when artificially inoculated.

All 24 samples subjected to paired-end RNA-seq revealed 1148.54 million clean reads with an average of 47.85 million clean reads in each sample (Supplementary Table S3). The fastq files of raw data are submitted at Sequence Read Archive database (SRA accession PRJNA616079). An average of 40.82 million reads were mapped on the Brassica napus genome, and on average, 80% of the reads were uniquely mapped in all samples. Considering FPKM thresh hold value of more than 1 in all three replications of resistant lines inoculated with PDA without fungus (R-mock) and susceptible lines inoculated with PDA without fungus (S-mock), 40,491 and 39,104 genes were expressed, respectively. Keeping the values ranging from 0.1 to 0.9 in every replication 2237 and 1891 genes were expressed in R-mock plants and S-mock plants, respectively, this low expression might be due to the stage-specific or tissue-specific nature of the genes. The Pearson correlation coefficient between the replications in both resistant lines and susceptible lines was calculated as high as 0.9 in most of the cases (Supplementary Figure 1), indicating the credibility of the RNA-seq method and samples. RNA-seq data were also validated with qPCR analysis with nine genes selected for qPCR, The Pearson correlation was calculated among FPKM values of RNA-seq data of same genes and qPCR relative expression for the resistant lines and was found to be higher as 0.99, assuring the reliability of RNA-seq analysis (Supplementary Table S4).

In resistant lines inoculated with disease (R) compared to R-mock lines, more genes were differentially expressed than in susceptible lines inoculated with disease (S) compared to S-mock lines. In total, 19,970 differentially expressed genes (DEGs) were found in R lines compared to R-mock lines (Figure 6A), and 3936 DEGs were found in S lines compared to S-mock lines at all time points (Figure 6B). Notably, after disease inoculation, there were a large number of DEGs in resistant lines compared to its mock line than in susceptible lines compared to its mock line. We also compared the DEGs between the R and S lines, and there were 14,735 genes in total at all three time points, including time-specific or DEGs sharing different time points in response to disease (Figure 6C). The large number of DEGs in R lines compared to S lines in response to disease inoculation indicates the tendency of resistant lines to cope with the disease.

Among RNA-seq for differentially expressed genes, we searched for the genes differentially expressed in our QTL regions, as these QTLs are already identified to be the regions containing genes affecting the disease resistance phenotype.

We compared R lines versus R-mock lines and R lines versus S lines, and there were more DEGs in the R vs. R-mock comparison compared to the R vs. S comparison. A total of 11,917 DEGs were there in the R vs. R-mock comparison, 6682 DEGs in the R vs. S comparison, and 8053 genes were expressed mutually in both comparisons at all time points (Figure 6D). We searched for the DEGs common in the R vs. R-mock and R vs. S comparisons at the same time point only (Figures 6E–G), and there were 550, 1172, and 2459 DEGs at 24, 48, and 96 hpi, respectively. In this study, we focused on the QTL positions identified in 2 years with the same confidence interval for DEGs. These QTLs were also identified in past studies as explained earlier. For the QTL region, there were 253 DEGs for SRA9a, 89 for SRC2a and 91 for SRC3a between R versus R-mock and 154 DEGs for SRA9a, 214 for SRC2a and 60 for SRC3a between R versus S at all three time points (Supplementary Table S5). Considering the genes only upregulated in one comparison, i.e., R vs. R-mock or R vs. S, may not contribute enough to the disease resistance trait. As the DEGs only in R vs. R-mock and not DEGs in comparison to S lines might not contribute enough to resistance and the DEGs only in R lines compared to S lines might be other varietal differences and not in response to disease. In common, 14 upregulated genes were found in SRA9a, 14 in SRC2a and 8 genes in SRC3a between both comparisons at different time points (Table 5). Most of the upregulated genes identified in QTL regions that were differentially expressed in both comparisons had been reported in previous studies to be involved in biotic or abiotic stresses. We included these 36 genes from three QTLs as putative candidate genes for SSR resistance.

The plant defense mechanism involves a wide range of responses in the recognition and handling of the disease (Wang et al., 2018), and a complicated network of genes was found to be expressed differentially in response to SSR within QTL regions. Different families of genes were upregulated in the R lines compared to the R mock and S lines.

One of the TIR-NBS-LRR class genes, BnaC03g05380D, present in the QTL region of SRC3a, was differentially expressed in R plants compared to R-mock plants at 48 hpi with a log₂ fold-change (log₂FC) of more than 1, and its expression was entirely absent in S plants at all time points with a log₂FC of almost 10 at each time point. Although there were many genes of this family present in our QTL region, only this gene was found with high relative expression in R plants. The only function of NBS-LRR proteins reported so far is disease resistance, whereas its role in resistance is yet be confirmed for most of the genes. The simplest model suggested for NBS-LRR protein function is that they act as receptors and bind the effector molecules secreted by pathogens. Proteolysis either specific or general is predicted to have its role in controlling the amplitude of defense response and the magnitude of cell death related to hypersensitive response (McHale et al., 2006). TIR-NBS-LRR class genes on other chromosomes were also predicted in other studies to be involved in SSR resistance in Brassica napus (Wu et al., 2016a, b).

Two zinc-induced facilitator 1 genes that are close together in the C3 chromosome of the SRC3a QTL were highly upregulated in the R lines compared to the R-mock and S lines. BnaC03g06010D was upregulated in all comparisons at all time points, whereas BnaC03g06020D was upregulated at 48 hpi in both comparisons. The difference in expression between the R and S lines was more obvious than the difference between the R and R-mock lines. Zinc-induced facilitator 1 belongs to one of the largest families of transporters present in plants, the Major Facilitator Superfamily (MFS) (Pao et al., 1998). The functional significance of a great number of MFS transporters is yet not determined in the fundamental physiological processes (Remy et al., 2013). An important function of MFS transporters Arabidopsis ZIFL1 was observed in regulating polar auxin transport and drought tolerance (Remy et al., 2013). ZIFL genes in wheat were also characterized for their expression in response to Zn and Fe deficiency (Sharma et al., 2019).

BnaC02g09450D, a nitrilase 4 gene, was also highly upregulated in R lines at all three time points within all comparisons present in the QTL SRC2a. With reference to their role in plant-pathogen interactions, nitriles are involved in hormone synthesis, nutrient assimilation and detoxification of exogenous and endogenous nitriles (Howden and Preston, 2009). In plants, nitriles are involved in defense pathways and in mostly cases are linked with cyanide metabolism (Howden and Preston, 2009). Nitrile compounds are also involved in the metabolism of cyanogenic glycosides and glucosinolates and both molecules are recognized defense molecules in plants that provide protection against pathogen attack (Howden and Preston, 2009).

BnaC02g48820D, an EXORDIUM-like 2 gene in QTL SRC2a, was highly upregulated in R lines at 48 and 96 hpi compared to S and R-mock lines. While screening for potential mediators of brassinosteroid the EXORDIUM gene was identified as upregulated (Coll-Garcia et al., 2004). Brassinosteroid metabolism is known to be activated differentially and adopt changes in plants when encounter abiotic stresses or bacterial, fungal and viral pathogens via activation of different mechanisms (Bajguz and Hayat, 2009).

BnaC02g00990D, an ATAF2 subfamily of NAC proteins, was highly upregulated in R lines at all three time points of disease inoculation, with a more obvious difference being observed in comparison to S lines than compared to R-mock. These proteins can form multiple protein complexes and have the ability to regulate a number of cellular processes during plant development stages and in response to stress (Puranik et al., 2012). The ATAF subfamily transcription factors act as the convergence point for biotic and abiotic stress signaling (Mauch-Mani and Flors, 2009). ATAF2 function was identified as a repressor for the proteins related to pathogenesis in Arabidopsis (Delessert et al., 2005).

Two UDP-glucosyl transferase genes were upregulated at 24 hpi in both comparisons. BnaA09g33440D, a UDP-glucosyl transferase 73C7, was more than 8 log₂FC at 24 hpi and was highly upregulated at 96 hpi. BnaA09g29790D, a UDP-glucosyl transferase 74B1, was more than 1 log₂FC compared to both susceptible and resistant mock lines at 24 hpi present in the SRA9a QTL. The overexpression of BnUGT74B1 in B. napus increased the aliphatic and indolic glucosinolates levels by 1.7 folds in leaves and in response to Sclerotinia infection showed less severe disease symptoms and tissue damage compared with the wild type control. UGT 73B3 and 73B5 responses to Pseudomonas syringae pv. tomato disease resistance have previously been reported and emphasize the importance of plant UGTs in plant-pathogen interactions (Langlois-Meurinne et al., 2005).

The other genes included in the QTL regions with upregulated expression (Table 5) in resistant lines at different time points are explained briefly in this study with their functions in disease resistance. The phylogenetic analysis in both distance-based and sequence-based analyses revealed an association of glycosyl hydrolase family 31 protein with defense enzymes (Opassiri et al., 2006). Aspartate is a precursor for the biosynthesis of several amino acids and derived metabolites that have unique and essential roles in plant development and defense response (de la Torre et al., 2014). Phosphoglycerate mutase-like family proteins are involved in glycolytic activity, which is essential for guard cell function in Arabidopsis (Zhao and Assmann, 2011). ARM repeat superfamily protein functions have been reported in protein sorting via ubiquitination, while the protein degradation plays an important role in achieving proteome plasticity in plants under environmental stress conditions to respond effectively (Sharma and Pandey, 2015). α/β-Hydrolase superfamily proteins are characterized to evolve complex and specialized chemical adaptations in response to varying biotic and abiotic conditions (Mindrebo et al., 2016). Plant S-adenosyl-L-methionine-dependent methyltransferases are the main enzymes in many metabolic pathways, including phenylpropanoid pathway (Joshi and Chiang, 1998), while reprogramming of the phenylpropanoid pathway is reported for SSR disease resistance in soybean (Ranjan et al., 2019). Gene repression via histone deacetylase plays a key regulatory role in plants to respond environmental stress (Luo et al., 2017). Polynucleotidyl transferase, a ribonuclease H-like superfamily protein, has an important role in mRNA deadenylation and mediation of stress responses (Walley et al., 2010). The ELMO/CED-12 family protein is required for the engulfment of dying cells and for cell migration (Gumienny et al., 2001). Many stresses are known to affect ribosomal protein gene transcript levels (McIntosh and Bonham-Smith, 2006). The roles of the KH-domain RNA-binding protein At5g53060 of Arabidopsis is reported in jasmonic acid signaling and biotic-induced stress responses (Thatcher et al., 2015). Invertases have a role in the establishment of plant defense responses as during the infection pathogens use plant sugars for their own needs and force plants to alter their sugar content and trigger their defense responses (Tauzin and Giardina, 2014). Pyruvate kinase is a catalyst in the irreversible synthesis of pyruvate and ATP, which are both involved in multiple biochemical pathways (Baud et al., 2007). ACC synthase is reported to be involved in ethylene signaling in roots and responsible for immediate inhibition of root elongation in response to pathogen-associated molecular patterns (Tsang et al., 2011). Amino acid transporters are required for resource allocation processes during plant growth and development (Ortiz-Lopez et al., 2000). DNA-binding proteins from tobacco plants are identified as a novel class of plant-specific regulatory factors that are involved in plant-virus interactions (Castelló et al., 2011). The plant tryptophan biosynthetic pathway is involved in the production of many secondary metabolites with diverse functions (Zhang et al., 2008).

We have also observed that flowering time has an effect on the disease resistance trait in our study and in some past studies (Wei et al., 2014; Wu et al., 2019), but in this study, we only investigated the genes directly involved in resistance. Under natural conditions, it may be predicted that early flowering plants shed their petals earlier on leaf stem attachment and provide a long moist period for disease to develop before ripening compared to late flowering plants. However, late flowering plants might be involved in disease escape under natural conditions because they do not allow petals to stay on them for an extended period and take less time for ripening due to increased temperature. However, in artificial inoculation, we inoculated both genotypes at the same time when termination of flowering was completed in both genotypes, which may reduce the effect of the environment on disease development and emphasize genetic factors that are involved in disease resistance. There may be a difference in developmental stage at the time of artificial inoculation that might have its role in disease resistance, as susceptible lines have completed their flowering slightly earlier than resistant lines and are ahead in developmental stage in the course of ripening. There may be a more vigorous network of genes involved in different pathways that were activated in the resistant lines during the time of inoculation that might help the resistant lines cope better with disease. The colocalization of the QTLs involved in flowering and SSR resistance was also reported, and the genes that accelerate flowering time in plants are closely linked to genes that increase susceptibility to the disease (Wu et al., 2019). Common genetic regions are involved in flowering time and SSR resistance was also reported in other past studies (Wei et al., 2014). The early flowering with increased disease susceptibility was observed in our study which is consistent with the previous findings (Wei et al., 2014; Wu et al., 2019). There is a strong need to identify the pathways or mechanisms affecting flowering and SSR resistance at the same time, and emphasizing the study for genetic linkage between these two traits and how these traits may be used for improvement in cultivar development with increased resistance. We have attempted to avoid the role of flowering in our identification of disease resistance genes by the selection of extreme lines for resistance and susceptibility from the same population with almost the same flowering time. This approach helped us to identify DEGs in the same growth period of both lines that might be involved directly in resistance to SSR.

Our study supports the narrative that SSR resistance is a quantitative trait controlled by the complex network of genes, as stated earlier (Wu et al., 2013). Acquiring the complete resistance genes or network of genes involved in resistance remains a challenge. Identifying a resistant source with their regulatory mechanism is a complex task, especially in polyploid crops. Although there was no common QTL identified in all 3 years in this study, the correlation among the phenotypic measurements was high in all 3 years and most of the QTLs identified in this study are in accordance with previous findings (Zhao et al., 2006; Wu et al., 2013, 2016a, 2019; Wei et al., 2014, 2016). However, the environmental fluctuations like temperature, humidity and other factors may alter the plants adaptation in response to disease, especially for such a disease with a very complex pattern of pathogen–plant interaction. It was noted that both lines selected for RNA-seq initiated flowering more than 2 weeks later in 2019 compared to 2017 (Figure 5A). Considering the common QTLs in two seasons and with other studies, we sought to identify the genes within these QTL regions that contribute to SSR resistance. Comparative RNA-seq is a reliable technique that we can use to cover the gene expression analysis of all the genes within QTL regions. With the help of DEGs, we were able to identify some of the genes that make little but significant contributions to SSR resistance. The genes identified in this study have shown their relation to biotic or abiotic stress in previous studies, highlighting the importance of these genes and supporting our findings. One of the identified gene in our present study, BnUGT74B1, was reported earlier in an overexpression study to be directly involved in increasing glucosinolates content and increased resistance against S. sclerotiorum infection (Zhang et al., 2015). This consistence with the previously reported gene suggested that the identified genes may play important roles for development of resistant cultivars by introducing the reported genes via overexpression or with a breeding strategy to enhance their expression. We suggest that these genes identified in our studies may have a direct role in resistance and need to be studied more with their specific role to counter the fungal disease S. sclerotiorum. The detailed study of all the DEGs in resistant and susceptible lines with additional experiments will be discussed in a future article to understand the network of genes in the whole genome responsible for disease resistance.