A total of 211 single seed descent (SSD)-derived recombinant inbred lines (RILs) were generated from the cross YS#24 × YS#58. The parental lines YS#24 and YS#58 are stable RILs selected from the Yangmai6 × Sonalika cross and advanced to the F₁₂ generation. Yangmai6 is a Chinese source of SB resistance, and Sonalika is a susceptible cultivar of Mexican origin that has been under cultivation in India for more than five decades/during the green revolution. The parents of the RILs used in this study are similar with respect to the agronomical and phenological traits but harbor different SB resistance QTLs (Kumar et al., 2009). To develop RILs, field trials in the crop season were conducted on the Agricultural Research Farm, BHU, Varanasi, whereas off-season nurseries were raised for generation advancement at Wellington, Tamil Nadu, India. The whole procedure of RIL development and its evaluation is mentioned in the flow chart (Supplementary Figure S1).

The 211 RILs (F₅ and F₆) were evaluated at three hot spots in India, namely, Agricultural Research Farm, BHU, Varanasi (25°18′N, 83°03′E); Borlaug Institute for South Asia, Samastipur, Pusa (25°57′N, 85°40′E), Bihar; and Uttar Banga Krishi Vishwavidyalaya, Coochbehar (26°19′N, 89°27′E), West Bengal, during two consecutive crop seasons 2013–14(F₅) and 2014–15(F₆). The F₇ generation was evaluated only at BHU in crop season 2015–16. At each location, all RILs were planted along with their parents (YS#24 and YS#58) in two replications following a randomized complete block design (RCBD). Each line was sown in two rows of 2 m, keeping row-to-row and plant-to-plant distances of 20 and 5 cm, respectively. To ensure disease build-up and spread, two rows of the susceptible genotype Sonalika was planted after every 20th row and in alleys along with the plots. The susceptible spreader rows served as an inoculum source for epidemic development in addition to the direct inoculation on the RILs (Saxesena et al., 2017). Planting at all sites was done between the 1st and 10th of December each year to coincide with the post-anthesis with higher temperatures, conducive to disease development (Chaurasia et al., 2000). The experimental plots were fertilized at 120 kg N₂, 60 kg P₂O₅, and 40 kg K₂O per hectare. A complete dose of K₂O and P₂O₅ was given at the time of sowing. The N₂ was given in three splits of 60, 30, and 30 kg per ha at sowing, 21 days after sowing (at first irrigation), and 45 days after sowing (at second irrigation), respectively. A total of five irrigations were given to maintain sufficient moisture in the field.

The highly aggressive B. sorokiniana isolate HDBHU (NABM MAT1; NCBIJN128877) was used to create an artificial epiphytotic. The pathogen was multiplied on sorghum grains, suspended in water (10⁴ spores per ml), and sprayed in the evening at the time of flag leaf emergence as described earlier (Joshi and Chand, 2002). Disease severity (DS) was scored using a double-digit scale, first at the beginning of anthesis (GS 63) and then at the end of anthesis (GS 69) and third at the late milk stage (GS 77) (Saari and Prescott, 1975). The percentage DS was calculated as D1/9 × D2/9 × 100, where D1 is the first digit referring to the vertical progress of disease and D2 is the second digit indicating the extent of leaf area affected. The AUDPC based on DS recorded at three growth stages was calculated following Shaner and Finney (1977). AUDPC  =   ∑ i = 1 n − 1 [ { ( Y i + Y i + 1 ) / 2 } × ( t ( i + 1 ) − t i ) ] , where Yi = disease level at time t_i, t _(i+1)−t_i = time (days) between two disease scores, n = number of dates when disease was recorded.

To increase the statistical power and reduce false positives, multiple bulks were selected independently from each of the three locations in both years (Zou et al., 2016). AUDPC was analyzed from the pooled data of disease severity in all the screened environments, and three resistant bulks with lower AUDPC (R-bulk1, R-bulk2, and R-bulk3) and three susceptible bulks with higher AUDPC (S-bulk1, S-bulk2, and S-bulk3) were prepared (Figure 1). Each group bulked for resistance and susceptibility is composed of 10 RILs.

Seeds of each RIL of the resistant and susceptible bulks (30R + 30S RILs) and their parents (YS#24 and YS#58) were grown separately in a greenhouse; seedling leaves were harvested 14 days after sowing for RNA isolation. Total RNA was extracted from the leaf tissues of each RIL using a standard TRIzol method (Catalog # 15596018) and treated with DNase I to remove residual DNA contaminants. RNA samples were purified using the RNeasy Mini Kit (Qiagen, Hilden, Germany) and the quantitative and qualitative estimation was performed using the NanoDrop 1000 (Thermo Fisher Scientific) and the Agilent Bioanalyzer 2100, respectively. An equimolar concentration of RNA from 10 resistant individuals of RILs was pooled together to make resistant bulk 1, and the same was done to prepare the other resistant and susceptible bulks. Before cDNA library preparation, we enriched the RNA samples for transcripts using the absolute mRNA Purification Kit (Agilent Technologies, Santa Clara, CA, United States). The cDNA libraries were constructed, and Illumina paired-end adapters and barcode sequences were ligated onto the cDNA fragments (Pootakham et al., 2014). The pooled libraries were sequenced at QTLomics Technologies (Bengaluru, India) using the Illumina Next Seq 500 platform (Illumina, San Diego, CA, United States) to generate 76 bp paired end (PE) sequence reads for parents and bulked samples (resistance and susceptible) in three biological replicates.

The data from Illumina Next Seq 500 was passed through Fast QC to check the quality of the reads (Babraham Bioinformatics - FastQC A, 2012). The ends of the reads with low quality, adaptor contamination, and low-quality regions were trimmed using the fastx-tool kit (http://hannonlab.cshl.edu/fastx_toolkit). Identification of SNPs was done; using a reference-based approach, the UniGene build63 was used as reference transcriptome (NCBI:ftp://ftp.ncbi.nhi.gov/repository/UniGene/Triticum_aestivum/Ta.seq.uniq.gz). Since UniGene represents only the genic portion, the complexity of the wheat genome was further reduced (Sidhu et al., 2015). The reads were aligned to transcript sequence build 63 using the Burrows–Wheeler Aligner (BWA) (Li and Durbin, 2009) using default parameters for PE libraries. Misaligned reads were removed from the unigene alignment and the unigenes with a high bulk frequency ratio (BFR) shortlisted. The allele frequency ratio of the resistant allele (allele of the resistant parent) in the resistant bulk compared with the resistant allele frequency in the susceptible bulk is known as the BFR. However, the genome locations of the SNPs were derived by alignment of shortlisted unigenes to the reference genome. We ran BLAST alignment and took the best hit throughout the transcript rather than the short read. The whole transcript alignment would give the best alignment to their respective genomes rather than the short read. The alignments were converted to one binary alignment/map (BAM) file per sample. On average, 80% of the RNA-Seq reads from each sample were able to align to the reference transcriptome (Pootakham et al., 2014). The SNPs for each sample were collected using the SAMtool sv0.1.18 as described by Li et al. (2009).

To identify SNPs linked to SB resistance, polymorphic markers between the parents were identified using a custom Perl script followed by BFR calculation for the bulks following Trick et al. (2012). The algorithm was implemented in Perl and R software. BFR values were calculated independently for the comparisons of three bulks (Bulk 1:S-bulk1, R-bulk1; Bulk 2: S-bulk2, R-bulk2; Bulk 3: S-bulk3, R-bulk3). BFR values and depth of each SNP for all the bulks and parents were calculated in the R program. Potential SNPs linked to the trait were selected based on SNPs with a BFR ratio of six as the minimum threshold in all three bulk replicates. The unigene that contain putative SNPs were aligned to a repeat masked reference genome (transcriptome) of wheat using BLASTN, and the best hit for each transcript was recorded. The BLAST result was segregated according to the chromosome number, and the chromosomal regions based on SNP density (twice or greater than the average SNP density) were further shortlisted. To identify SNPs that were enriched for the corresponding parental allele, BFR was calculated in the appropriate bulk, YS #24 derived SNPs for the resistant bulks, and YS#58 derived SNPs for the susceptible. Finally, to classify and prioritize the SNPs across each bulk, the frequency of the allele (SNP index) calculated at each SNP position (Takagi et al., 2013) was estimated, and then the ratio between the bulks (BFR) for each SNP was determined (Trick et al., 2012). Thus, a high BFR in the resistant parent (YS#24) derived SNP was indicative of an allele that is very frequent in the resistant bulk while depleted in the susceptible. In addition, a threshold of BFR >6 was set to select putative SNPs for the presence or absence of polymorphism between bulks for further validation as done previously in wheat (Ramirez-Gonzalez et al., 2015).

The tetra-primer amplification refractory mutation system-PCR (ARMS-PCR) is a simple and economical technique that produces an allele-specific reaction (Ye et al., 2001; Ruiz-Sanz et al., 2007) The principle behind this technique is that the allele-specific primer amplifies a region specific to the base present at the 3′terminus, thus making it allele-specific (Newton et al., 1989; Ye et al., 2001). The genomic region flanking 300 bp on both ends of putative SNPs was extracted and formatted using the Perl script, and the SNPs with the highest score were selected to design primers. The batch primer 3 web tools (http://probes.pw.usda.gov/batchprimer3/) were finally used to create the tetra-primer ARMS-PCR (You et al., 2008).

The genomic DNA of RILs used for BSR-Seq was extracted from the leaves of young seedlings using a cetyltrimethylammonium bromide (CTAB) protocol (Saghai-Maroof et al., 1984). The DNA pellet was vacuum dried, dissolved in DNase-free water, and stored at −20°C. A target-specific tetra-primer ARMS PCR amplification of all the resistant and susceptible bulks along with the parents was performed. The PCR was performed in a total volume of 10 μl reaction mixture, containing 30 ng DNA, 1X PCR buffer [75 mM Tris-HCl (pH9.0), 50 mM KCl, 20 mM (NH₄)₂SO₄], 0.33 pmol of each outer primer (forward and reverse), 0.5 pmol of inner forward primer, 0.83 pmol of inner reverse primer, 2.5 mM MgCl₂, 0.125 mM of each dNTPs, 1U of Taq Polymerase (3B DNA polymerase, 3B Black Bio Biotech India Ltd.). PCR reactions were performed in an Agilent SureCycler 8800 using the following program: initial denaturation at 95°C for 5 min, 38 cycles consisting of denaturation at 95°C for 1 min, annealing at 62°C for 90 s, extension at 72°C for 2 min, followed by a final extension at 72°C for 10 min. The PCR product was resolved on 3.5% agarose gel and stained with ethidium bromide to visualize.

To perform the functional analysis of the transcripts, Blast2GO v 2.5 was used (Conesa et al., 2005). It is a Gene Ontology–based annotation tool and was found effective in the functional characterization of sequence data (Conesa and Götz, 2008). For functional characterization of the transcripts, we performed NCBI-BLASTX analysis using transcripts homologous to annotated proteins in the nr database with the criterion of E-value (threshold of 1e-03) and alignment size (threshold length 33). Furthermore, the transcript sequences were categorized for gene ontology (i.e., GO terms) into three groups: molecular function, biological process, and cellular component. The pathways for the selected transcripts were also delineated using the Kyoto encyclopaedia of genes and genomes (KEGG) database.

In all environments, the AUDPC of the resistant parent YS#24 was consistently and considerably higher than that of the susceptible parent YS#58 (Table 1). The mean AUDPC values of the resistant (YS#24) and susceptible parents (YS#58) were 336.29 and 820.55, respectively, over the environment (Table 1). Based on the pooled data, the RIL population for SB AUDPC revealed considerable phenotypic variation, ranging from 231.70 to 836.80. The RILs also exhibited a high coefficient of variation (CV) for AUDPC across environments, ranging from 8.80% (BHU14) to 16.92% (UBKV14), whereas the broad-sense heritability (h ²) values of AUDPC were lowest at BISA14 (0.65) and highest at BHU14 (0.87). Therefore, resistant and susceptible bulks were constituted by taking the samples independently using extreme phenotypes (resistant and susceptible) from both tails as illustrated in Figure 1.

The bulk samples utilized for BSR-Seq produced 429.40 million raw reads, containing 32,634.40 million base pairs. The number of raw reads in the sequenced samples ranged from 30.40 million (7.08%) in the resistant bulk-3 to 171.70 million (39.99%) in the susceptible parent (Table 2). The highest mapping of reads to the reference genome was obtained for the susceptible parent (171.70 million reads), followed by susceptible bulk-3 (43.90 million reads), resistant parent (42.50 million reads), resistant bulk-1 (37.90 million reads), and susceptible bulk-2 (36.50 million reads), and the lowest values were obtained for resistant bulk-3 (30.40 million reads) (Figure 2A). All the raw sequence reads have been deposited at NCBI, and the provisional Sequence Reads Archive (SRA) identifiers were obtained under the accession codes SRR5948908.

Each sample yielded an average of 80% high-quality (free of adaptor contamination and low-quality areas) RNA-Seq reads that were matched to the wheat reference transcriptome. Susceptible bulk-1 had the highest percentage alignment (80.38% with 535390 SNPs), followed by susceptible bulk-3 (80.09% with 613480 SNPs), and resistant bulk-1 had the lowest percentage alignment (79.01% with 669279 SNPs) (Table 2). The resistant parent had the highest number of SNPs (843836), whereas the minimum number (535390) was in the susceptible bulk-1 with the highest percentage of alignment ( Figures 2B,C ).

A total of 1,379,122 SNPs were identified on 91,855 transcripts of wheat samples used in the present investigation (Table 3). The putative SNPs linked to SB were selected based on SNPs with BFR >6 in all three bulk samples, i.e., 3666 SNPs present on 1837 transcripts in Bulk-1 (S-bulk1: R-bulk1), 3635 SNPs on 2056 transcripts in Bulk-2 (S-bulk2: R-bulk2), and 2130 SNPs on 1443 transcripts in Bulk-3 (S-bulk3: R-bulk3). A total of 7860 SNPs with >6 BFR were found across all bulks (Bulk-1: Bulk-2: Bulk-3) on 3950 transcripts. Out of 7860 SNPs, only 1055 SNPs were present on 506 transcripts detected as homozygous and polymorphic between parents (Table 3). The transcripts carrying homozygous and polymorphic SNPs with other details are listed in Supplementary Table S1.

The distribution of trait-linked SNP markers in the A, B, and D genomes for each homoeologous group and percentage of the total number of SNPs is given in Table 4. The 1055 polymorphic SNP markers, bar plotted on 21 wheat chromosomes, formed seven unevenly distributed homoeologous groups (Figure 3). The highest number of SNP markers were identified in the B genome (532 SNP; 50.42%), followed by A (311 SNP; 29.46%) and D (212 SNP; 20.09%). The number of SNPs per linkage group ranged from 16 (2D) to 198 (5B) (Figure 4). The homologous group 5 had the highest markers (313 SNP, 29.67%), followed by group 3 (207 SNP, 19.60%), and group 4 had the lowest (79 SNP, 7.49%) (Table 4).

The maximum number of polymorphic SNP with BFR >6 was found on chromosome 5B (198 SNP), followed by 3B (136 SNP) (Figure 4). Thus, a total of 334 SNPs of 3B and 5B chromosomes were present on 142 transcripts, out of which 60 SNPs were found only on five transcripts. Out of these five, one transcript of the 3B chromosome had eight SNPs, and four transcripts of the 5B chromosome had 52 SNPs (Table 5). Among the four transcripts of 5B, a maximum of 27 SNPs were present on transcript gnl/UG/Ta#S61812294, followed by transcripts gnl/UG/Ta#S17985740 (19 SNPs) and gnl/UG/Ta#S61830716 (4 SNPs), while transcript gnl/UG/Ta#S65715070 had the lowest number (2 SNPs) (Table 5). The chromosomal distribution of 60 polymorphic SNP present on five different transcripts of the 3B and 5B chromosomes associated with SB resistance is shown in Figure 5, which indicates their relative position.

In this study, allele-specific tetra-primer ARMS PCR primers were designed to develop an assay for SB resistance. The SNPs on transcript gnl/UG/Ta#S61799095 of 3B chromosomes and transcript gnl/UG/Ta#S61830716 and gnl/UG/Ta#S17985740 of chromosome 5B were used to design the tetra-primers for ARMS PCR (Supplementary Table S2). A primer (Ta_S61830716_1262_3) from chromosome 5B amplified 142 bp fragments in both the resistant parent and resistant bulk-1, whereas a 142 bp fragment was not amplified in the susceptible parent, but in a few samples of susceptible bulk-1, a faint band appeared in lane 10–12 (Figures 6A,B). The primer (Ta_S61830716_1262_3) showed clear differentiation between resistant and susceptible genotypes.

The transcript of 3B chromosomes (gnl/UG/Ta#S61799095) shared homology with acetyl-CoA acetyltransferase, whereas two transcripts of 5B (gnl/UG/Ta#S17985740 and gnl/UG/Ta#S61830716) shared homology with proteinase/protease and the remaining two with phospholipase C1 and exohydrolase proteins (Supplementary Table S3). A total of 346 GO terms were identified for the selected transcripts and further categorized into molecular function, biological process, and cellular component (Supplementary Figure S2). The maximum number of GO terms fall into the molecular function category (231 terms), which were further grouped into hydrolase activity, transferase activity, peptidase activity, and catalytic activity. There were 70 terms associated with the biological process, which fall into the category of protein catabolic process, carbohydrate metabolic process, metabolic process, phospholipid catabolic process, and fatty acid beta-oxidation (Supplementary Figure S2). A few terms were associated with cellular components, viz., peroxisome, membrane, an integral component of the membrane, lysosome, and extracellular space. Furthermore, as shown in Supplementary Figure S3, the most highly identified pathways in the study were fatty acid metabolism, valine, leucine, isoleucine metabolism, and benzoate degradation.