The spot-blotch-resistant wheat genotype IC566637 was obtained from the Indian Council of Agricultural Research (ICAR)—National Bureau of Plant Genetic Resources (NBPGR), New Delhi, India, while Agra Local served as the susceptible control. Seeds were sown in 4-inch. pots containing pre-sterilized soil and maintained at 25 ± 2 °C in the greenhouse at the Division of Plant Pathology, ICAR—Indian Agricultural Research Institute (IARI), New Delhi. A pure culture of B. sorokiniana isolate BS-22 (virulent on Agra Local, NCBI accession number ON892001, Indian Type Culture Collection identification number 9114), mass-cultured using potato dextrose agar medium and autoclaved sorghum grains (25–30 days), was used for the study. Seedlings at the three- to four-leaf stage were inoculated with the spore suspension (10⁴ conidia/mL) using a precision hand atomizer. Sterile distilled water was used as untreated control. To maintain high humidity (>90%), the plants were kept in a humidity chamber box for 24–48 h. Treatments consisted of susceptible inoculated (SI) and susceptible untreated control (SC) as well as resistant inoculated (RI) and resistant untreated control (RC). Leaf samples from inoculated and control plants were collected at 0, 12, 24, 36, 48, and 60 h post-inoculation in liquid nitrogen and stored in a deep freezer at -80°C until further use. Each sample consisted of three biological replicates.

For small RNA sequencing (sRNA-seq), two biological replicates per treatment were used, with samples pooled across time points. Total RNA was extracted from leaf tissue using the Direct-ZOL mini kit (Zymo Research) following standard protocols. RNA concentration and purity were measured using a Nanodrop spectrophotometer (Thermo Scientific, 2000). The proportion of miRNA in the samples were determined using a Bioanalyzer (Agilent), while RNA quantification was performed with the Qubit RNA HS assay kit (Thermo Scientific). QIAseq^® miRNA Library Kit (Qiagen) standard protocol was used to generate the sRNA-seq libraries. In brief, 63 ng of total RNA was used as initial input. The 3′ adapters were selectively ligated to the 3′OH group of miRNAs, followed by 5′ adapter ligation. The ligated fragments were subjected to reverse transcription with Unique Molecular Index (UMI) assignment by priming using reverse transcription primers. The resulting cDNA was then amplified and barcoded through PCR amplification (17 cycles). The Illumina-compatible sequencing libraries were quantified using the Qubit fluorometer (Thermo Fisher Scientific, MA, USA). Finally, single-end sequencing was performed for 50 cycles on Illumina NovaSeq 6000 High Output sequencing platform using SE50 read length and sequencing chemistry.

Raw reads of 16–40 bases in length were filtered and mapped to the Triticum aestivum reference genome (PRJNA669381) using Bowtie (https://bowtie-bio.sourceforge.net/index.shtml) (Langmead et al., 2009). The reads were then specifically aligned to sRNA. The obtained reads were further aligned with noncoding RNA database to filter out rRNA, tRNA, and snoRNAs. The remaining unaligned and high-quality reads (Q30 ≥ 90%) were classified as known or novel miRNAs through a homology approach against Viridiplantae miRNAs from miRBase22 (miRBase) (Friedländer et al., 2012) using the NCBI-BLAST-2.2.304 with an e-value cutoff of e^-4 and non-gapped alignment. Non-homologous sequences were analyzed using Mireap_0.22b, and those forming stem loop structures were considered as novel miRNAs. To eliminate the potential effect of B. sorokiniana miRNAs on wheat, all of the detected miRNAs were cross-referenced against B. sorokiniana genome.

Read counts for known and novel miRNAs were obtained to analyze their expression patterns. DE analysis was conducted using DESeq (version 2) (Anders and Huber, 2010), an R package with a threshold of log₂ (fold change, FC) ≥1 and false discovery rate (FDR) ≤0.05. Read variations were normalized using DESeq2’s library normalization method, where size factors were calculated and each count was divided by its corresponding factor. The average normalized read counts of samples in each condition were utilized for differential gene expression calculation. For target prediction analysis, miRNAs with a copy number ≥5 were selected. These miRNA sequences were used as input along with wheat–B. sorokiniana transcriptome sequences (unpublished) to the miRanda tool (Griffiths-Jones et al., 2008), with the following default parameters: (1) using the strict alignment in the seed region (offset positions 2–8), preventing detection of target sites containing gaps or non-canonical base pairing in this region, (2) maximum expectation, and (3) target accessibility—allowed maximum energy to unpair the target site (UPE) 25. The DE miRNAs across eight sRNA-seq samples were clustered based on their log₂FC values, with those displaying similar expression patterns grouped using R package Heatmap version 1.0. The associated target genes of DE miRNAs were analyzed for the GO and the KEGG enrichment to identify significantly enriched functional categories.

The identified target gene sequences were annotated with the BLAST2GO software (BLAST: Basic Local Alignment Search Tool (nih.gov), GO (http://www.geneontology.org/), KEGG (http://www.genome.jp/kegg/), Pfam (http://pfam.xfam.org/), InterPro (InterPro (ebi.ac.uk)), and Panther (pantherdb.org). Furthermore, the identified target genes were analyzed for GO (http://www.geneontology.org/) and KEGG (http://www.genome.jp/kegg/) (Kanehisa and Goto, 2000) enrichment analysis with FDR cutoff of 0.05. The networks of DE miRNAs along with their associated DE target genes were prepared using Cytoscape (v.3.10.2) (Shannon et al., 2003).

The same samples used for sRNA-seq were also employed to validate miRNA expression and their corresponding target genes. RNA isolated from infected leaves was converted into cDNA using the Verso cDNA synthesis kit (Thermo Scientific) with miRNA-specific stem loop primers (Kramer, 2015). As per the protocol described by Varkonyi-Gasic et al. (2007), forward and universal reverse primers were designed for the selected nine miRNAs. The U6 small nuclear RNA (U6 snRNA) was the internal control gene for normalization. The correlation of miRNA expression in qRT-PCR (log₂FC) and sRNA-seq was analyzed using the scatter plot. Based on the sequences of nine predicted target genes associated with the nine selected miRNAs, their specific qPCR primers were designed on GeneScript (Real Time PCR Primer Design - Real Time PCR Probe Design - GenScript). β-actin was used as an internal control gene. The qRT-PCR analysis of miRNAs and their corresponding target genes was performed using the Bio-Rad CFX96 system (Bio-Rad Laboratories, Inc., India) in the genomics laboratory of the Discovery Centre, ICAR-IARI, New Delhi, India. The list of primers used is given in Table 1.

An average of 19.73 million reads were generated as raw reads from the eight wheat samples comprising untreated and B. sorokiniana-inoculated wheat leaves. This data has been uploaded in the Sequence Read Archive (SRA) database of NCBI with Bioproject number PRJNA944096. After stringent quality filtering and removal of low-quality reads, adapters, and sequences shorter than 16 nt, an average of 12.2 million high-quality reads (Q30 ≥ 90%) with a read length between 16 and 40 nt were retained for downstream analysis (Supplementary Table S1). A total of 736 miRNAs were identified across all samples, including 418 known and 318 novel miRNAs (Supplementary Table S2). The majority of the identified miRNAs were 21 nt (41.57%) in length, followed by 22 nt (21.60%), 20 nt (13.72%), and 24 nt (11.27%) (Figure 1). The symptomatic response of the resistant genotype IC566637 and suscpetible variety Agra local at 4 days post-inoculation is shown in Figure 2.

Among the 736 identified miRNAs across the control and treated groups of the two wheat lines, 140 were DE miRNAs in response to infection. The number of DE miRNAs was slightly higher in SI than in RI (Supplementary Table S3). Specifically, 44 miRNAs were upregulated and 43 downregulated in SI, whereas 37 were upregulated and 47 downregulated in RI. A total of 31 DE miRNAs were common between RI and SI, indicating that some miRNA-mediated regulatory responses are common to both resistant and susceptible backgrounds. However, genotype-specific expression patterns were also evident: 30 and 23 miRNAs were uniquely up- and downregulated in RI, while 25 and 31 were unique to SI, respectively. A total of 12 miRNAs were found to be upregulated in both the RI and SI, with all exhibiting higher log₂FC values in SI, suggesting stronger transcriptional activation of these miRNAs under susceptible background (Figures 3A, B). Interestingly, several miRNAs showed contrasting expression in RI and SI; tae-miR1122b-3p, tae-miR9662a-3p, bdi-miR5054, gma-miR395i, osa-miR5072, and ppt-miR894 were downregulated in RI and were found to be upregulated in SI (Figure 4A). Conversely, ata-miR5181-3p, osa-miR5503, cme-miR156j, osa-miR5076, and tae-miR7757-5p were found to be upregulated in RI and downregulated in SI (Figure 4B).

For the 140 DE miRNAs, a total of 894 putative target genes were predicted and functionally annotated based on GO, KEGG, Panther, InterPro, and Pfam (Supplementary Table S4). The identified target genes represented a broad range of functional categories. Majority of the identified target genes were associated with cellular signaling processes such as protein kinases, including serine/threonine protein kinases and calcium-transporting ATPases. Several targets were involved in plant-defense-associated functions, including xylanase inhibitors (TAXI-IV and XI-IB), beta-1,3-glucanase (Glb3), peroxidases (Prx111-A and Prx112-M), chitinases, and resistance-related proteins such as RPP13-like, Nucleotide-binding Apaf-1 R proteins CED-4 (NB-ARC), and LRR proteins. In addition, a number of target genes encoded transcription factors, including WRKY10 and zinc finger proteins (ZFPs), which are often linked with regulatory control of gene expression. Several transporter proteins were also identified among the targets, such as members of the sulfite exporter TauE/SafE family, putative zinc transporter (ZIP1), and ABC transporters. Furthermore, multiple transferase enzymes were predicted, including acyl transferases (acT1) and ubiquitin transferases (RING-type E3).

The GO enrichment analysis of predicted target genes was grouped into biological process, cellular components, and molecular functions (Figure 5). For the specific DE miRNAs in IC566637, the predicted target genes participated in 144, 114, and 12 GO terms for biological process, molecular function, and cellular component, respectively (Supplementary Table S5; Figure 5A). However, the predicted target genes for DE miRNAs specific in Agra Local participated in 92, 104, and five GO terms for biological process, molecular function, and cellular component, respectively (Supplementary Table S6; Figure 5B). The KEGG pathway enrichment analysis indicated that the predicted target genes were mainly involved in 15 significantly enriched pathways (Table 2). In the resistant genotype IC566637, the most enriched pathways included biosynthesis of secondary metabolites, valine, leucine, and isoleucine biosynthesis, diterpenoid biosynthesis, fatty acid metabolism, and the MAPK signaling pathway. In the susceptible genotype Agra Local, the enriched pathways included linoleic acid metabolism, phenylpropanoid biosynthesis, glucosinolate biosynthesis, cysteine and methionine metabolism, valine, leucine, and isoleucine metabolism, plant–pathogen interaction, and ubiquinone and other terpenoid–quinone biosynthesis.

Hierarchical clustering and heat map analysis revealed distinct expression patterns with upregulated (green) and downregulated (red) miRNAs in IC566637 and Agra Local, forming two separate clades that indicate expression heterogeneity (Figure 6). For IC566637, clade A contained upregulated and clade B with downregulated miRNAs in RI compared to RC samples (Figure 6A). However, in Agra Local, clade A shows downregulated and clade B with highly upregulated miRNAs in SI compared to SC samples (Figure 6B). The predicted DE target genes of miRNA for IC566637 in clade A were predominantly enriched in GO terms with transporter activity (calcium ion transmembrane transport, ABC transporters, ATPase coupled transporter activity, transmembrane transporters, and methyltransferase), cytokinin and hormone metabolic process, RNAi effector complex, RISC complex, and COPII vesicle coat and vesicle tethering complex (Figure 6C). However, the predicted DE target genes in IC566637 were mainly enriched with organic acid, carbohydrate, carboxylic acid metabolic process, extracellular region, lyase and acyltransferase activity, hydrolase activity acting on glycosyl bonds, and sulfur compound binding (Figure 6D). The predicted DE target genes of miRNA in clade A of Agra Local were mainly enriched in carbohydrate metabolic process, ion transport, small molecule metabolic process, organic substance transport, cell wall, external encapsulating structure, extracellular region, and transporter activity (Figure 6E). For miRNAs in Agra Local (clade B), the DE-predicted target genes were mainly enriched in hormone-mediated signaling pathway, cellular response to hormone stimulus, endogenous stimulus and organic substance, pre-ribosome small subunit precursor, eukaryotic translation elongation factor 1 complex, nuclear membrane, CCAAT-binding factor complex, SLIK (SAGA-like) complex, clathrin-coated vesicle, 1-phosphatidylinositol binding, minor groove of adenine–thymine-rich DNA binding, DNA-binding transcription activator activity, and transcription corepressor activity (Figure 6F). The results suggest that miRNAs exhibiting varying expression levels could have unique functions in regulating their target genes engaged in different biological pathways during B. sorokiniana and wheat interaction.

In this study, we evaluated the expression of nine miRNAs selected based on lowest p-value and highly variable expression (either significantly upregulated or downregulated) and their corresponding target genes associated with defense response against B. sorokiniana in IC566637 and Agra Local. Among these, the expression of four known miRNAs (tae-miR9662a-3p, osa-miR5072, ata-miR9863b-3p, and ata-miR9863b-5p) and five novel miRNAs (xxx-m0030-3p, xxx-m0005-3p, xxx-m0007-5p, xxx-m0009-5p, and xxx-m0031-5p), respectively, was assessed at five different (12, 24, 36, 48, and 60 h) time points using qRT-PCR (Figure 7). The expression of these miRNAs varied considerably across different time points. Notably, three known miRNAs (tae-miR9662a-3p, osa-miR5072, and ata-miR9863b-3p) exhibited contrasting expression patterns: they were downregulated in IC566637 and upregulated in Agra Local. For instance, the expression of tae-miR9662a-3p reduced considerably in RI, which corresponded with the upregulation of its target gene (XM_044510133.1: gene involved in defense response to fungus) that reached 8.06-fold increase at 36 h. Conversely, in SI, the expression of tae-miR9662a-3p drastically increased to 59.74-fold at 48 h, while its target gene was downregulated. Similarly, osa-miR5072 and ata-miR9863b-3p were found to be downregulated in RI and upregulated in SI. The osa-miR5072 exhibited its highest expression in SI at 29.12 times that of the control, while ata-miR9863b-3p reached 3.59 times that of the control in SI. Correspondingly, the target genes for osa-miR5072 (NB-ARC protein) and ata-miR9863b-3p (NB-LRR protein) were downregulated in SI while being upregulated in RI, reaching up to 17.89-fold at 24 h and 2.83-fold at 36 h, respectively. The expression of ata-miR9863b-5p was significantly upregulated in both RI and SI, reaching up to 4.11-fold and 6.07-fold at 48 h, respectively. Notably, its target gene, putative disease resistance (RPP13-like protein 1), showed differential regulation. In RI, the gene was positively regulated by the miRNA with an upregulation of 13.19-fold compared to the control. Conversely, in SI, a negative regulation was observed, where the target gene was found to be downregulated. Among the predicted novel miRNAs evaluated for expression, four miRNAs—xxx-m0030-3p, xxx-m0005-3p, xxx-m0007-5p, and xxx-m0009-5p—were found to be downregulated in SI. Contrarily, the corresponding target genes for xxx-m0030-3p (cyclin-dependent protein serine/threonine kinase), xxx-m0005-3p (peroxidase), and xxx-m0007-5p (ATP-binding, involved in defense response) were up-regulated in SI, rising to 5.23-fold at 60 h, 3.3-fold at 60 h, and 1.56-fold at 24 h, respectively, showing a negative regulatory relationship. However, for the target gene for xxx-m0009-5p (protein kinase), a positive regulation was observed, where the target gene was also downregulated in SI. The other miRNA xxx-m0031-5p was significantly upregulated in SI, rising to 5.12-fold compared to the control at 48 h, with corresponding downregulation of its target gene (transcription cofactor). The expression profile of these DE miRNAs in qRT-PCR showed a strong correlation (R² = 0.915) with sRNA-seq data (Figure 8). The stem loop structures of validated miRNAs were depicted (Figure 9).