Brachypodium distachyon genotype Bd21-3 was selected for this study. All plants, including control plants, were cultured in a growth chamber at 22°C under a 16/8 h photoperiod. Bd21-3 is compatible with M. oryzae race RO1-1. Approximately 5 weeks after sowing, the plants were inoculated with fresh M. oryzae spores (1 × 10⁵ spores mL^–1) by spraying (Wang et al., 2020). Mock inoculation was performed in parallel with sterile water. Inoculated plants were placed in a moist chamber at 28°C in complete darkness for the first 24 h. Inoculated leaves and mock-treated control leaves were sampled at 0, 24, and 48 h post-inoculation (hpi) and snap-frozen in liquid nitrogen. Three independent biological replicates were used for each treatment.

Total RNA was extracted from B. distachyon leaves using TRIzol (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s instructions. RNA concentration was determined using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, MA, United States). Approximately 1 μg of total RNA was used to construct a sRNA library following the protocol of the mRNA-Seq sample preparation kit (Illumina, San Diego, CA, United States). Single-end sequencing (50 bp) was performed using an Illumina Hiseq2500 at the LC-BIO (Hangzhou, China) (Zhang et al., 2019).

Raw reads were processed with ACGT101-miR (LC Sciences, Houston, TX, United States) to remove adapters, junks, repeats, low-complexity sequences, and common non-coding RNA families (rRNA, tRNA, snRNA, and snoRNA). The remaining clean reads with lengths of 18--24 nts were blasted to miRBase 21.0¹ to identify known miRNAs. Novel miRNA candidates were blasted to the genome of B. distachyon. The pre-miRNAs were mapped to specific species genomes to determine their genomic locations. The unmapped sequences were analyzed using BLAST against the B. distachyon genome. Then, hairpin RNA structures were predicted from the flank 120 nt sequences using RNAfold² as previously described (Meyers et al., 2008). All raw sequencing data were deposited into the NCBI Short Read Archive under the BioProject: PRJNA751253 (Accession number: SRR15316624–SRR15316628).

Approximately 20 μg of total RNA was used to prepare the degradome library (Zhai et al., 2020). Poly(A) + RNA was bound to mRNA capture beads as input RNA. Fifty adapters were ligated to RNAs containing 50-monophosphates, reverse transcribed, and PCR amplified. The purified cDNA library was used on Illumina’s Cluster Station and sequenced on an Illumina Hiseq 2500 at the LC-BIO. Raw sequencing reads were used to analyze potentially cleaved targets using the CleaveLand 3.0 pipeline.³ Differentially expressed miRNAs revealed by the normalized deep-sequencing read counts were analyzed by Student’s t-test. The threshold for significance was set at p < 0.05 and for high significance at p < 0.01 in each test.

The sequences of B. distachyon were obtained from Ensembl Plants⁴ genome-centric portal for plant species. We assigned genes to different categories based on the predicted functions of their Arabidopsis and rice homologs. BLAST2GO (Conesa et al., 2005) was then used to obtain Gene Ontology (GO) annotations at default settings of Fisher’s exact test (p < 0.05), false discovery rate (FDR) correction using the Hochberg method, and five minimum numbers of mapping entries against a species-specific pre-computed background reference. The KEGG13 pathway was analyzed using the ClueGO plug-in14 and Cytoscape software V2.8.215 to identify significant pathways of the differentially expressed genes (DEGs) (Bindea et al., 2009).

Inoculated and control leaves were harvested at 0, 16, 48, and 120 hpi. First-strand cDNA was synthesized using a quantitative real-time PCR (qRT-PCR) system (Takara, Shiga, Japan) according to the manufacturer’s instructions. All qRT-PCR was performed on a CFX96 Real-Time System (Bio-Rad, Munich, Germany) using SYBR Green I (Thermo Fisher Scientific, Carlsbad, CA, United States). The total volume of PCR reaction was 25 μL. PCR conditions were as follows: 94°C for 5 min, 35 cycles (94°C for 30 s, 56°C for 30 s, and 72°C for 30 s) and 72°C for 5 min. UBC18 was used as an internal control for qRT-PCR analysis to normalize the data (Hong et al., 2008). The relative expression of R-genes was quantified using the comparative 2^–ΔΔCT method (Voegele and Schmid, 2011). Primers for defense genes were designed using the DNAMAN 9 software (Supplementary Table 1). Three independent biological replicates were used for each treatment.

The region containing the cleavage site of BdNAC21 was cloned into the pBI121 vector containing the β-glucuronidase (GUS) reporter gene. The precursor of miR164c was integrated into the pBI121 vector to replace the GUS gene. Strain GV3101 of Agrobacterium tumefaciens was cultured in LB medium at 28°C to introduce the recombinant vectors into Nicotiana benthamiana leaves. Spectra of A. tumefaciens strains were normalized to an equivalent optical density of 0.8 at 600 nm (OD600) before injection. The strains mixed with pBI121-BdNAC21 and pBI121-miR164c were infiltrated into N. benthamiana leaves (Wang et al., 2017). Each treatment was infiltrated with 1 mL of the liquid containing the strains. Histochemical staining and quantitative GUS detection were performed was performed using 4-methylumbelliferyl-β-D-glucuronide (4-MUG) as the substrate (Jefferson et al., 1987).

β-Glucuronidase activities were measured using a Gemini XSP spectrophotometer plate reader (Molecular Devices, Sunnyvale, CA, United States) at wavelengths of 365 and 455 nm. Each assay was then performed in duplicate.

The phenotype of B. distachyon infected with M. oryzae (isolate RO1-1) was identified. Large disease lesions were observed in susceptible B. distachyon plants (Bd21-3) 5 days after inoculation. The early stage in the infection of B. distachyon with M. oryzae is reportedly the key time point of the plant–fungus interaction. To validate whether sRNAs mediate the response of B. distachyon to M. oryzae infection, we examined the dynamic changes in miRNA expression in B. distachyon infected by M. oryzae at 0, 24, and 48 hpi using Illumina Hiseq2500. Fifteen sRNA libraries were constructed including M-0 h, M-24 h, I-24 h, M-48 h, I-48 h (mock and inoculation is abbreviated as M and I respectively; M-0 h: Mock 0 hpi; M-24 h: Mock 24 hpi; I-24 h: Inoculation 24 hpi; M-48 h: Mock 48 hpi; I-48 h: Inoculation 48 hpi; three biological replicates were performed per treatment). A total of approximately 13.6 million reads were obtained from each library (from 10,265,735 to 15,611,967 reads per library) (Figure 1A and Supplementary Table 2). After removing low-quality reads, unique reads (from 1,564,288 to 2,568,854 reads per library) were obtained (Figure 1B). The length distribution of the sRNAs ranging from 18 to 25 nt was further analyzed. sRNAs 24 nt in length were the most frequent, followed by sRNAs 21 nt in length (Figure 2A). We found that four types of nts had an average total miRNA nt bias (Figure 2B). Moreover, the first position in the 18–25 nt miRNAs was predominantly uracil (U), except for 24 nt miRNAs starting with adenine (A) (Figure 2C). Unique reads were annotated against miRBase and Rfam to eliminate rRNA, tRNA, and snoRNA (Supplementary Table 2).

A total of 432 known miRNAs were identified in the sRNA libraries by analyzing the unique clean reads that were aligned with miRBase (Figure 3A). These known miRNAs belonged to 77 miRNA families (Supplementary Table 3). We found that the miR156, miR164, miR166, and miR167 families contain more than 30 members. These known miRNAs were further compared to the data of 12 different plant species including Arabidopsis thaliana and Oryza sativa to explore their evolutionary characteristics (Figure 3B). As shown in Figure 3C, the 77 known miRNA families exhibited different numbers of homologous sequences in most of the comparisons. Fifteen known miRNAs were found to be conserved (Figure 3C). The remaining sequences that did not match any of the known miRNAs but mapped to the genome were then further analyzed based on whether they formed hairpins. A total of 288 predicted candidate miRNAs (PC miRNAs) were identified (Supplementary Table 4). Venn diagrams were created based on the comparison of five library miRNA sequences. The results showed that 26 and 12 miRNAs were only found in M-0 h and M-24 h groups, respectively. Moreover, seven miRNAs were unique to the I-48 h group (Figure 4).

Differentially Expressed miRNAs (DEMs) examined at 24 hpi (M-24 h: Mock 24 hpi; I-24 h: Inoculation 24 hpi) and 48 hpi (M-48 h: Mock 48 hpi; I-48 h: Inoculation 48 hpi), yielded the absolute value of log2 fold change higher than 1.5 and a p-value below 0.05. We found seven DEMs at 24 hpi; six miRNAs were upregulated, whereas one miRNA was down-regulated (Figure 5A). Compared with M-48 h, we identified 19 DEMs in I-48 h, including four upregulated miRNAs and 15 downregulated miRNAs (Figure 5B and Supplementary Table 5). The expression of eight novel miRNAs changed at 48 hpi. Compared with M-0 h, we found 25 and 28 DEMs in B. distachyon infected by M. oryzae at 24 and 48 hpi, respectively (Supplementary Figure 1). Moreover, the number of DEMs between I-24 h and I-48 h was 25 (Supplementary Figure 1).

To identify the target genes of miRNAs involved in the B. distachyon–M. oryzae interaction, a mixed degradome library was generated from the M-0 h, M-24 h, M-48 h, I-24 h, and I-48 h groups. Approximately 99.59% (24,200,004) sequenced raw reads of the degradome library were mappable when compared with genomic DNA; the number of unique mappable reads was 5,944,614. A total of 19,431,096 (79.96%) unique reads were mapped to the transcripts of B. distachyon protein-coding genes (Supplementary Table 6). Furthermore, strict conditions (category ≤2 and raw reads ≥10) were applied to screen the target genes to achieve accurate results for miRNAs and target genes in the degradome data library. A total of 2,126 genes were identified as targets for 308 miRNAs, including 258 known miRNAs and 49 novel miRNAs. The conserved miRNAs have multiple targets including transcription factors. Among the most conserved miRNAs and targets, we detected miR156a-R1-squamosa promoter-binding-like protein (SPL), miR164-NAC, miR396-growth-regulating factor (GRF), and miR444-MADS-box transcription factor (MADS). Novel miRNAs and targets were also detected, such as PC-3p-201628_13-THO2 (Supplementary Table 7).

To better understand the regulatory roles of DEMs and their targets in the interaction between B. distachyon and M. oryzae, GO enrichment annotations were performed at 24 and 48 hpi. GO enrichment analysis revealed that there were significantly enriched processes involved in plant defense responses, such as JA and ethylene-dependent systemic resistance. Furthermore, acyl-CoA oxidase and kinesin-binding activities were among the significantly enriched processes at 24 hpi (p < 0.05) (Figure 6A). At 48 hpi, ubiquitin-protein transferase activity and kinase activity were among the significantly enriched processes (Figure 6B). Moreover, to elucidate the exact biological process by which DEGs might participate in plant–fungus interactions, we performed a KEGG enrichment analysis. Our results revealed that “peroxisome” and “ubiquitin mediated proteolysis” were significantly enriched pathways involved in plant defense responses at 24 and 48 hpi, respectively (Supplementary Figure 2 and Supplementary Table 8).

To validate the miRNA/target regulatory unit, the changes in expression of miRNAs and mRNAs during the B. distachyon–M. oryzae interaction were analyzed using qRT-PCR. Five miRNAs, including one novel miRNA and four known miRNAs (Figure 7) were randomly selected to show cleavage events via degradome sequencing. DEMs and their corresponding targets were selected for qRT-PCR analysis. Consistent with the analysis of sequencing data for miRNAs, five target genes were negatively regulated by the corresponding miRNAs. Moreover, miR164c/Bd4g02060 (BdNAC21), miR169-p5/Bd1g21177 (BdMED12), and miR9494/Bd5g06390 (BdFKBP65) were mainly regulated at 24 hpi. However, miR5182-p3/Bd3g41480 (BdPDH2) and PC-3p-15479/Bd3g28280 (BdPRP19) were more significantly regulated at 48 hpi than at 24 hpi.

To confirm the suppression of BdNAC21 by miR164c, we performed transient co-transformation of N. benthamiana (Figure 8A). The recombinant vector pBI121 contained GUS as a reporter gene. Leaves infiltrated with pBI121-BdNAC21 showed the GUS staining. In leaves infiltrated with mixture of pBI121-BdNAC21 and pBI121-miR164c, GUS staining was reduced. The amount of GUS in each leaf was also assessed to further validate results of the histochemical observations (Figure 8B). We found that the fluorescence from the empty vector and pBI121-BdNAC21 was enhanced on extending the reaction time for the inoculated leaf samples. No fluorescence was observed in leaves injected with pBI121-miR164c. However, the fluorescence of the leaf injected with the mixture strain (pBI121-BdNAC21 and pBI121-miR164c) exhibited a slow increase. These results further confirm that BdNAC21 is the gene targeted by miR164c.