Comparative Analysis of Transcriptome and sRNAs Expression Patterns in the Brachypodium distachyon—Magnaporthe oryzae Pathosystems

The hemibiotrophic fungus Magnaporthe oryzae (Mo) is the causative agent of rice blast and can infect aerial and root tissues of a variety of Poaceae, including the model Brachypodium distachyon (Bd). To gain insight in gene regulation processes occurring at early disease stages, we comparatively analyzed fungal and plant mRNA and sRNA expression in leaves and roots. A total of 310 Mo genes were detected consistently and differentially expressed in both leaves and roots. Contrary to Mo, only minor overlaps were observed in plant differentially expressed genes (DEGs), with 233 Bd-DEGs in infected leaves at 2 days post inoculation (DPI), compared to 4978 at 4 DPI, and 138 in infected roots. sRNA sequencing revealed a broad spectrum of Mo-sRNAs that accumulated in infected tissues, including candidates predicted to target Bd mRNAs. Conversely, we identified a subset of potential Bd-sRNAs directed against fungal cell wall components, virulence genes and transcription factors. We also show a requirement of operable RNAi genes from the DICER-like (DCL) and ARGONAUTE (AGO) families for fungal virulence. Overall, our work elucidates the extensive reprogramming of transcriptomes and sRNAs in both plant host (Bd) and fungal pathogen (Mo), further corroborating the critical role played by sRNA species in the establishment of the interaction and its outcome.


Introduction
The ascomycete fungus Magnaporthe oryzae (Mo) (anamorph: Pyricularia grisea) is the causal agent of rice blast, one of the most devastating and widespread diseases of cultivated rice, reducing yields up to 30% annually [1][2][3]. Members of the Magnaporthe genus can also infect a variety of other cereals, including barley, rye and wheat, making Mo a major threat to global food security [4,5]. Currently, blast management strategies rely on a combination of fungicide applications (e.g., azoles), development of resistant cultivars, and agronomic practices such as removal of crop residues, water management and crop/land rotation [6,7].
In the early stage of foliar infection Mo behaves as a biotroph, forming a biotrophic interfacial complex (BIC) between the primary invasive hyphae (derived from the penetration peg) and the infected host cell, where it secretes small molecules (effectors) to modulate the interaction [8,9]. Subsequently the fungus forms secondary hyphae and spreads to neighboring cells, undertaking a lifestyle change switching to a necrotrophic growth, with the appearance of the characteristic blast lesions on leaves [5]. Mo is able to infect all aerial parts of rice including nodes, panicles and necks, and has been shown to produce necrotic lesions on both rice and barley roots, although its lifestyle in roots seems to depend on the inoculation method [4,10]. Mo infections have also been established The Magnaporthe oryzae (Mo) 70-15 genome encodes three putative AGO and two DCL proteins, previously identified by domain search and phylogeny with known Neurospora crassa (Nc) RNAi machinery genes [25][26][27]. The corresponding protein sequences were obtained from NCBI: XP_003714515.1 (MGG_01541T0, MoDCL1), XP_003715365.1 (MGG_12357T0, MoDCL2), XP_003716704.1 (MGG_14873T0, MoAGO1), XP_003717504.1 (MGG_13617T0, MoAGO2) and XP_003714217.1 (MGG_01294T0, MoAGO3); they were included in new phylogenetic trees to corroborate previous findings ( Figure S1A,B). To assess the conservation of key AGO and DCL domains, prediction was carried out with Simple Modular Architecture Research Tool (SMART). All three AGOs have conserved domain structures with five characteristic domains: N-terminal domain, DUF1785, PAZ (Piwi Argonaut and Zwille), L2 and PIWI (P-element Induced WImpy testis) ( Figure S1C; Table S1). Additionally, multiple sequence alignment (MSA) confirmed the conservation of the DEDD catalytic tetrad (Asp/Glu/Asp/Asp) and the QF-V motif (Glu/Phe/Val) in the MoAGO PIWI domains ( Figure S2). Likewise, both MoDCLs have the four conserved domains required for the cleavage of dsRNAs: DEXDc, HELICc, dicer-dimer and RIBOc ( Figure S1D; Table S1). Subcellular localization of MoAGOs and MoDCLs was assessed using PSI-predictor; MoDCL1, MoAGO1 and MoAGO3 were predicted to primarily localize in the nucleus, while MoDCL2 and MoAGO2 were predicted to primarily localize in the cytosol. MoAGO1 and MoAGO2 also had secondary predicted localizations, the first in the cytosol and the second in plastids (Table S2).
Protein interaction analysis with STRING did not highlight any experimentally proven interaction for MoAGO and MoDCL proteins, but probable interactions could be inferred from data derived from homologs, found to either interact or co-express with several proteins in other species (Table S3A,B). In particular, all MoAGOs were predicted to interact with the two MoDCL proteins, a U5 small nuclear ribonucleoprotein component (MGG_13500) and a cell cycle control protein cwf14 (MGG_06309). For MoAGO1, additional interaction was predicted with Pumilio-family RNA-binding repeat protein (MGG_03158), ATP-dependent RNA helicase DED1 (MGG_02762), high-affinity glucose transporter (MGG_13651) and four uncharacterized proteins. In addition to interaction with AGO proteins, both fungal MoDCLs were predicted to interact with ATP-dependent DNA helicase MPH1 (MGG_04429), 30S ribosomal protein S16 (MGG_02598), WD domaincontaining protein (MGG_06727) and three uncharacterized proteins.
3D protein structure modeling for both MoDCL and MoAGO was performed with SWISS-MODEL. While no acceptable QMEAN Z-scores (>−4) values were obtained for the two DCLs and AGO2 models due to the lack of suitable reference structures, two models for AGO1 (model 1 based on hAGO1, GMQE = 0.52 QMEAN = −3.98 and model 3, based on hAGO2, GMQE = 0.51 QMEAN = −3.31) and model 6 for AGO3 (based on hAGO4 template, GMQE = 0.46 QMEAN = −3.50) passed the first quality check and were further validated with PROCHECK and WHATCHECK. Model 1 of AGO1 had a higher percentage of residues in the core region of the Ramachandran plot (88.4%) and a better Ramachandran Z-score (−1.808) compared to model 3 (87.6% and −1.531 Z-score), while the AGO3 model scored 86.7% and −2.088, respectively. As a result, model 1 of AGO1 and model 6 of AGO3 were selected as the best models and visualized with PyMOL ( Figure 1A,B). validated with PROCHECK and WHATCHECK. Model 1 of AGO1 had a higher percentage of residues in the core region of the Ramachandran plot (88.4%) and a better Ramachandran Z-score (−1.808) compared to model 3 (87.6% and −1.531 Z-score), while the AGO3 model scored 86.7% and −2.088, respectively. As a result, model 1 of AGO1 and model 6 of AGO3 were selected as the best models and visualized with PyMOL ( Figure  1A,B). Next, we assessed the virulence of Δmoago1, Δmoago2, Δmodcl1, Δmodcl2 and the double mutant Δmodcl1/2 using Bd21-3 as a host. Negative effects on vitality of the five mutants vs. Mo wildtype (wt) could not be detected when inspecting conidia germination and appressorium formation on poly-L-lysine coated glass slides by 24 h after the onset of germination ( Figure S3). In contrast, virulence of these mutants was compromised both on Bd leaves and roots ( Figure S4). Since we expected differences in infection phenotypes depending on test conditions, three different leaf inoculation protocols were tested: i.) upon drop inoculation of detached leaves all mutants showed significant (Kruskal-Wallis test, p < 0.05) reduction in necrotic lesion size, ranging from a −20% to −48% reduction of relative necrotic area compared to Mo wt (Δmoago1 −48%; Δmoago2 −32%; Δmodcl1 −39%; Δmodcl2 −24%; Δmodcl1/2 −20%; Figure 2A); ii.) upon spray inoculation of detached leaves, virulence of Δmoago1 and Δmodcl1 was virtually unaffected as they showed the typical necrotic lesions at 5 days post inoculation (DPI), compared to Δmoago2, Δmodcl2 and Δmodcl1/2, showing significantly reduced necrotic lesion areas (compared to Mo wt: Δmoago2 −40%; Δmodcl2 −83%; Δmodcl1/2 −68% necrotic lesion sizes; Figure 2B); iii.) similar phenotypes were observed in the whole seedling spray assay, with the exception of Δmodcl2 and Δmodcl1/2 that appeared avirulent and did not show any detectable lesions (reduction in the necrotic lesion size compared to Mo wt: Δmoago2 −79%; Δmodcl2 −100%; Δmodcl1/2 −100%; Figure 2C). Finally, mutants were inoculated on Bd roots and fungal presence was measured at 5 DPI by qPCR based on the fungal actin. All mutants showed a significant reduction in fungal biomass, with Δmoago2, Δmodcl2 and Δmodcl1/2 being the most severely affected (reduction in MoActin levels compared to wt: Δmoago1 −70%; Δmoago2 −78%; Δmodcl1 −47%; Δmodcl2 −87%; Δmodcl1/2 −87%; Figure 2D). These results confirm and extend previous findings that Mo's RNAi machinery is required for disease development in cereals and, in particular, that DCL2 is a key factor not only for sRNA biogenesis but also for fungal virulence. Next, we assessed the virulence of ∆moago1, ∆moago2, ∆modcl1, ∆modcl2 and the double mutant ∆modcl1/2 using Bd21-3 as a host. Negative effects on vitality of the five mutants vs. Mo wildtype (wt) could not be detected when inspecting conidia germination and appressorium formation on poly-L-lysine coated glass slides by 24 h after the onset of germination ( Figure S3). In contrast, virulence of these mutants was compromised both on Bd leaves and roots ( Figure S4). Since we expected differences in infection phenotypes depending on test conditions, three different leaf inoculation protocols were tested: (i.) upon drop inoculation of detached leaves all mutants showed significant (Kruskal-Wallis test, p < 0.05) reduction in necrotic lesion size, ranging from a −20% to −48% reduction of relative necrotic area compared to Mo wt (∆moago1 −48%; ∆moago2 −32%; ∆modcl1 −39%; ∆modcl2 −24%; ∆modcl1/2 −20%; Figure 2A); (ii.) upon spray inoculation of detached leaves, virulence of ∆moago1 and ∆modcl1 was virtually unaffected as they showed the typical necrotic lesions at 5 days post inoculation (DPI), compared to ∆moago2, ∆modcl2 and ∆modcl1/2, showing significantly reduced necrotic lesion areas (compared to Mo wt: ∆moago2 −40%; ∆modcl2 −83%; ∆modcl1/2 −68% necrotic lesion sizes; Figure 2B); (iii.) similar phenotypes were observed in the whole seedling spray assay, with the exception of ∆modcl2 and ∆modcl1/2 that appeared avirulent and did not show any detectable lesions (reduction in the necrotic lesion size compared to Mo wt: ∆moago2 −79%; ∆modcl2 −100%; ∆modcl1/2 −100%; Figure 2C). Finally, mutants were inoculated on Bd roots and fungal presence was measured at 5 DPI by qPCR based on the fungal actin. All mutants showed a significant reduction in fungal biomass, with ∆moago2, ∆modcl2 and ∆modcl1/2 being the most severely affected (reduction in MoActin levels compared to wt: ∆moago1 −70%; ∆moago2 −78%; ∆modcl1 −47%; ∆modcl2 −87%; ∆modcl1/2 −87%; Figure 2D). These results confirm and extend previous findings that Mo's RNAi machinery is required for disease development in cereals and, in particular, that DCL2 is a key factor not only for sRNA biogenesis but also for fungal virulence. Detached second youngest leaves of three-week-old Bd seedlings were drop inoculated with 10 µL Mo conidia suspension 50 × 10 3 spores mL −1 in 0.002% Tween20 and kept under high humidity at 16 h light/8 h dark cycle at 22 • C/18 • C. Fungal pathogenicity was assayed via ImageJ software at 5 DPI. The experiment was conducted two times (n = 8 plants per experimental group) with similar results. Comparisons between groups was performed via ANOVA and Tukey's range test for multiple comparisons. (B) Detached Bd leaves were spray inoculated with a total of 250 µL Mo conidia suspension 50 × 10 3 spores mL −1 in 0.002% Tween20 and kept and evaluated as in (A). The experiment was conducted two times (n = 8 plants per experimental group) with similar results. Comparisons between groups was performed via Welch one-way test coupled with pairwise t-tests with Benjamini-Hochberg p-value adjustment. (C) Three-week-old Bd seedlings were spray inoculated with Mo conidia suspension 120 × 10 3 spores mL −1 in 0.002% Tween20 and kept and evaluated as in (A). The experiment was conducted two times (n = 18 plants per experimental group) with similar results. Comparisons between groups was performed via Kruskal-Wallis test and Dunn's test of multiple comparisons. (D) Roots of seven-day-old Bd seedlings were inoculated with 1 mL of Mo conidia suspension 125 × 10 3 spores mL −1 in 0.002% Tween20 for 3 h. Afterwards, seedlings were transplanted in small pots (3 cm Ø) and grown for an additional 5 days before harvesting. Fungal amount was calculated by qPCR based on the ratio of fungal actin (MoActin). The experiment was conducted two times (n = 6 roots per experimental group) with similar results. Comparisons between groups was performed via ANOVA and Tukey's range test for multiple comparisons. (A-D) Letters represent statistical difference among all groups' means (α < 0.05). Asterisks represent statistical difference of each group against wildtype (wt) (* p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001).
Minor overlaps were observed comparing downregulated Bd genes between the biological setups, with only two shared between the foliar infections, and five between the leaf (4 DPI) and the root setup ( Figure 5A). Among the upregulated genes, the highest overlap was seen between the foliar infections, sharing 134 DEGs, compared to 16 DEGs shared between root and leaf 2 DPI, and 72 between root and leaf 4 DPI ( Figure 5B). Interestingly, 15 Bd genes were upregulated in all setups, including defense-related ABC transporter (BdiBd21-3.3G0465100.1), anthranilate synthase component II (BdiBd21-3.5G0159100.1), protein kinase xa21 (BdiBd21-3.3G0144800.1) and secologanin synthase-like (BdiBd21-3.2G0563800.1) ( Table 2).   Figure S5). A notable overlap of fungal DEGs was detected between the two foliar infections, with 346 up-and 498 downregulated genes shared, compared to 220 up-and 172 downregulated DEGs shared between leaf 2 DPI and root infection, and 498 up-and 203 downregulated DEGs shared between leaf 4 DPI and root infection. Overall, 142 genes were significantly downregulated in all setups, while 168 were consistently upregulated in all samples, including genes involved in fungal development, metabolism and pathogenicity ( Figure 5C,D).
Expression of DCL and AGO genes was assessed in both organisms in leaves and roots, which showed that BdDCL1 and MoDCL2 had the highest expression among the DCL family, and BdAGO3 and MoAGO3 among the AGO family ( Figure 6). samples, including genes involved in fungal development, metabolism and pathogenicity ( Figure 5C,D).
Expression of DCL and AGO genes was assessed in both organisms in leaves and roots, which showed that BdDCL1 and MoDCL2 had the highest expression among the DCL family, and BdAGO3 and MoAGO3 among the AGO family ( Figure 6).
Interestingly, BdAGO9, which is the closest homolog to AtAGO1 [35], was detected expressed in all three setups and had the highest read count of the AtAGO1-like clade. Moreover, BdAGO1, BdAGO5 and BdAGO15 were only expressed in the root samples, while BdAGO4, BdAGO7, BdAGO11, BdDCL3a and BdDCL4 were not detected in any of the datasets. Finally, the majority of DCL and AGO genes were not differentially regulated during the infection, with the exception of BdDCL3b and BdAGO3 (minor downregulation; >−1 log2FC) and BdDCL1-2a and MoAGO3 (minor upregulation; <1 log2FC) (Table S4).

Gene Ontology Enrichment (GOE) and Defense-Related Gene Expression in Mo-Infected Bd
GOE analysis using AgriGO v2 was performed on all Bd DEGs identified in leaves and roots by mRNAseq. In all three datasets, GO:0003824 (catalytic activity), GO:0016491 (oxidoreductase activity) and GO:0044710 (single-organism metabolic process) were enriched (Table S5). The highest number of Gene Ontology (GO) terms enrichment-in particular those related to various cellular metabolic processes-was recorded in leaves at 4 DPI, which corresponds to the necrotrophic fungal growth phase ( Figure S6). Consistent with the extensive metabolic reprogramming highlighted in the GOE analysis, the highest number of DEGs also was observed in these samples. Specifically, we detected significant upregulation of known Bd defense response genes such as pathogenesis-related 1 (PR1), PR3 and PR5, receptor-like protein kinases and Myb and WRKY transcription factors ( Table 2). Only minor overlap was observed between defense-related genes expression patterns in the leaf 2 DPI and root samples. Three predicted PR proteins were upregulated in both 4 DPI and root setups, including PR10 and PRB1-2 (log2FC > 3 in both, while not induced at Interestingly, BdAGO9, which is the closest homolog to AtAGO1 [35], was detected expressed in all three setups and had the highest read count of the AtAGO1-like clade. Moreover, BdAGO1, BdAGO5 and BdAGO15 were only expressed in the root samples, while BdAGO4, BdAGO7, BdAGO11, BdDCL3a and BdDCL4 were not detected in any of the datasets. Finally, the majority of DCL and AGO genes were not differentially regulated during the infection, with the exception of BdDCL3b and BdAGO3 (minor downregulation; >−1 log 2 FC) and BdDCL1-2a and MoAGO3 (minor upregulation; <1 log 2 FC) (Table S4).

Gene Ontology Enrichment (GOE) and Defense-Related Gene Expression in Mo-Infected Bd
GOE analysis using AgriGO v2 was performed on all Bd DEGs identified in leaves and roots by mRNAseq. In all three datasets, GO:0003824 (catalytic activity), GO:0016491 (oxidoreductase activity) and GO:0044710 (single-organism metabolic process) were enriched (Table S5). The highest number of Gene Ontology (GO) terms enrichment-in particular those related to various cellular metabolic processes-was recorded in leaves at 4 DPI, which corresponds to the necrotrophic fungal growth phase ( Figure S6). Consistent with the extensive metabolic reprogramming highlighted in the GOE analysis, the highest number of DEGs also was observed in these samples. Specifically, we detected significant upregulation of known Bd defense response genes such as pathogenesis-related 1 (PR1), PR3 and PR5, receptor-like protein kinases and Myb and WRKY transcription factors (Table 2). Only minor overlap was observed between defense-related genes expression patterns in the leaf 2 DPI and root samples. Three predicted PR proteins were upregulated in both 4 DPI and root setups, including PR10 and PRB1-2 (log 2 FC > 3 in both, while not induced at 2 DPI) and two were found significantly upregulated in the two foliar datasets, including PR1 (with a log 2 FC > 6 in these datasets and not induced in the roots).

GOE and Gene Expression in Mo during Bd Infection
Next, significantly up-and downregulated Mo genes were subjected to GOE analysis, showing the highest number of enriched GO terms in the root samples, including several terms relating to the interaction with the host (Table S6; Figure S7). Overall, many infectionrelated genes were significantly upregulated, including effectors AVR Pita1 and PWL2 (Pathogenicity toward Weeping Lovegrass) and NOXs (superoxide-generating NADPH oxidases) ( Table 3).

sRNA Reprogramming in Bd and Mo at Early Infection Stages
We assessed the expression of sRNAs in the same biological material utilized for mRNAseq. Before high-throughput next generation sequencing (NGS), multiplexed sRNA libraries were size selected for 15 to 35 nt fragments (140-160 nt including TruSeq adapters).
Single end sequencing on Illumina HiSeq1500 platform generated between 22 million (mil) and 38 mil reads each (Table S7). Quality check of raw reads was performed with FastQC, adapters were removed with cutadapt and reads were size selected between 18 to 32 nt. The organism of origin of the trimmed reads was predicted by mapping via Bowtie alignments to both Bd21-3 and Mo 70-15 genomes [44] (Bd21-3 v1.1 DOE-JGI, http://phytozome.jgi.doe.gov/). Ambiguous reads that could not be assigned to the organism of origin with high confidence were excluded to avoid miscalling. As expected, most reads in Mo-infected samples were assigned to Bd (with 100% match) and not to the fungus (with at least two nucleotide mismatches) (Table S7). Size distribution of genome matched unique sRNA reads followed a similar trend throughout samples, with the Mo reads showing a minor peak around 20 nt and Bd reads a considerable peak at 24 nt ( Figure 7A,B). Fungal unique sRNA reads were then compared among datasets derived from Mo axenic culture and Mo-infected plant tissues and classified as either exclusive or shared between samples ( Figure 8A). Some 5708 Mo-sRNAs were identified in Mo-infected roots, of which 3263 (57.15%) were found exclusively in the infected sample and not in the axenic culture. A total of 7215 Mo-sRNAs were identified in infected leaves during the biotrophic phase (2 DPI), of which 4399 (60.97%) were found exclusively in infected samples. Finally, a total of 63,017 were found in infected leaf tissue during the necrotrophic phase (4 DPI), of which 46,212 (73.33%) were exclusively found in infected samples. Samples for sRNA sequencing by Illumina HiSeq1500 were taken from different setups: leaf biotrophic phase ("2 DPI"), leaf necrotrophic phase ("4 DPI") and root ("root"). . Samples for sRNA sequencing by Illumina HiSeq1500 were taken from different setups: leaf biotrophic phase ("2 DPI"), leaf necrotrophic phase ("4 DPI") and root ("root").
Next, we compared unique Bd-sRNA reads in leaf and root samples from Mo-infected and mock-treated Bd ( Figure 8B). A large number of sRNAs were found in infected tissues: 571,644 in leaves at 2 DPI, of which 326,657 (72.24%) were specific for infection; 415,469 at 4 DPI, of which 265,172 (69.06%) were specific for infection; and 597,158 Bd-sRNAs in infected roots, of which 346,259 (77.92%) were specific for infection. These data show that the majority of unique sRNAs of both interacting organisms are expressed exclusively during the interaction, suggesting that they are relevant to the establishment and outcome of the disease. We further selected unique Bd-sRNAs that (i.) were either found exclusively in infected plant tissues or (ii.) showed higher reads in infected vs. mock-inoculated tissue. Interestingly, the size distribution of these induced Bd-sRNA did not show a change in length preference compared to the total unique reads ( Figure 7C,D). Relative size distribution (in percentage) of unique filtered sRNA reads assigned to (A) Mo or (B) Bd. Reads were assigned to either Mo or Bd only if aligning 100% to the organism of origin genome and had at least two mismatches to the interacting organism genome. (C,D) Relative size distribution of unique filtered sRNA reads assigned to (C) Mo or (D) Bd and induced or increased in infected samples compared to controls (i.e., axenic fungal cultures and non-inoculated plants, respectively). Samples for sRNA sequencing by Illumina HiSeq1500 were taken from different setups: leaf biotrophic phase ("2 DPI"), leaf necrotrophic phase ("4 DPI") and root ("root").

Infection-Related Bd miRNAs
Filtered Bd-sRNA reads were analyzed with Shortstack 3.8.5 to identify infectionrelated upregulation of potential plant miRNAs. The program identified 14 putative miRNA-generating loci in leaves at 2 DPI, 15 at 4 DPI and 16 in the roots. Matching predicted miRNA generating clusters in the uninfected samples were identified via their genomic coordinates, and numbers of total reads per cluster were compared between infected and uninfected samples to select exclusively expressed or upregulated miRNA loci. The selected clusters' sequences were aligned to known Bd miRNA stem loop precursor sequences obtained from the miRBase database: three known miRNAs were upregulated in leaves at 2 DPI, six at 4 DPI and four in the roots ( Table 4). Structures of the miRNA precursors are shown in Figure S8. miRNA-generating clusters were identified with Shortstack in all datasets. Genomic coordinates were used to detect and compare the same loci in infected and control samples, and only clusters exclusive or increased in the infected samples were selected. Clusters and mature miRNAs were compared to known miRNA/miRNA precursor from miRBase. Abbreviation: RPM = reads per million.

In Silico Prediction of Bd-sRNAs Targeting Mo Transcripts
Given that plant-derived sRNAs can move into fungal pathogens during fungal infections [39,40], we further explored the possibility that Bd-sRNAs have corresponding downregulated Mo transcripts. In total, 20-22 nt Bd-sRNAs ( Figure 9) that (i.) originate from the non-coding regions of the Bd21-3 genome, and (ii.) show a higher read count in the Mo-infected compared to non-infected sample were selected. We found 424, 302 and 681 Bd-sRNA in leaves at 2 DPI and 4 DPI and in roots, respectively, with predicted targets in the Mo transcriptome (Table 5). In the next step, we selected for targets found downregulated in the mRNAseq analysis, which further reduced the predicted Mo ck-sRNAs to 314 in leaves at 2 DPI, 236 at 4 DPI and 263 in roots, respectively, corresponding to potentially targeted mRNAs downregulated in leaves at 2 DPI (484), 4 DPI (527) and roots (183) ( Table 5). Downregulated Mo targets include cell wall-related genes such as chitin deacetylase 1 (MGG_05023), cell wall protein MGG_09460 and virulence genes, such as CAP20 (MGG_11916) ( Table 7). By comparing predicted Mo targets in infected tissues, we found a considerable overlap between leaf and root samples (42 Mo mRNAs) and between the two leaf samples (140 Mo mRNAs) ( Table 7; Figure 10B).  We searched the Pathogen-Host Interactions Database (PHI-base) for available information on loss-of-virulence mutants for the predicted Mo target genes. A list of downregulated predicted Mo targets and corresponding PHI-base phenotypes is shown in Table 8. Of note, the predicted Mo genes are involved in virulence and pathogenicity, including CON7 transcription factor (MGG_05287) predicted to be targeted in the Moinfected root, the effector AvrPiz-t (MGG_09055) predicted to be targeted in the leaf at 4 DPI and the vacuolar glucoamylase SGA1 (MGG_01096) predicted to be targeted in the leaf at 2 DPI. Additionally, potential targets shared between more than one setup included (MGG_00620), Sso1 (MGG_04090) encoding a SNARE protein and YHM2 (MGG_07201) encoding mitochondrial DNA replication protein (2 DPI leaf vs. root), as well as MoATG17 (MGG_07667) encoding an autophagy-related protein (biotrophic vs. necrotrophic leaves), whose respective mutants are also known to be compromised in pathogenicity [45,46] (Tables 7 and 8).

Discussion
Based on the work done to elucidate the function and 3D structure of Brachypodium distachyon AGO and DCL proteins published by Šečić et al. [35], we set out to investigate these RNAi machinery components in Magnaporthe oryzae. Mo has being extensively studied as a model for RNAi in fungi, and it is known to encode three AGOs, two DCLs and three RdRPs [25,26]. Interestingly, while the phylogeny with other fungal RNAi machineries such as N. crassa, Mucor circinelloides, Cryphonectria parasitica and Schizosaccharomyces pombe and the domain conservation have been reported, we found no information regarding the proteins' predicted interactome, subcellular localization and 3D protein structures, which would help elucidate the function and relevance of these proteins in vivo. We confirmed that all three fungal AGOs have the characteristic PAZ (Piwi Argonaut and Zwille) and PIWI (P-element Induced WImpy testis) domains: the first is involved in the recognition of the 3 end of the guide sRNA molecule [47], while the second is an RNase H domain that confers AGOs the ability to cleave single-stranded RNAs complementary to the guide sRNA sequence [48]. A closer look at the PIWI domains via MSA confirmed the conservation of both the QF-V motif, involved in the sorting of sRNAs based on their sequence and secondary structure [49] and the DEDD catalytic tetrad, critical for the protein 'slicer' activity [48]. Similar to NcDCLs, both MoDCLs lacked predicted PAZ domain but differed from the N. crassa ortholog showing predicted dicer-dimer domains and lacking predicted dsRNA-specific ribonuclease and ds-RNA binding (DSRM) domains [27]. At this point of the analysis, no major differences were recorded among either Mo protein families, which would help elucidate the differences in both sRNA production and virulence reported in the RNAi KO mutants by Raman et al. [50]. Interestingly, protein subcellular localization helped shed additional light on MoAGO3 and MoDCL2 non-redundant roles within their families, with MoAGO3 predicted to localize exclusively in the nucleus, compared to MoAGO1, predicted in both the nucleus and the cytosol, and MoAGO2, predicted to localize in the cytosol and plastid. Moreover, we predicted MoDCL2 to exclusively localize in the cytosol, while MoDCL1 in the nucleus, confirming the distinct roles for these proteins in contrast with their redundant function N. crassa [51]. Additionally, predicted interactomes of MoAGOs showed significant differences between MoAGO1 vs. MoAGO2 and MoAGO3, which would point to a unique role of this protein. Specifically, the predicted interaction with an ATP-dependent RNA helicase DED1 (MGG_02762) and a Pumilio-family RNA-binding repeat protein (MGG_03158) would suggest a unique role for MoAGO1 in translation regulation [52,53]. KO mutants of these RNAi components have been previously characterized both for their sRNA accumulation patterns and their virulence on barley [27,50]. Except for ∆moago3, which did not sporulate and could not be further utilized, we tested the morphology and virulence of the available mutants in our newly established Bd pathosystems. While all five RNAi mutants showed unaltered growth, morphology and conidiation in axenic cultures, we detected significant differences in their ability to infect Bd tissues. Interaction between plants and its host pathogen is a dynamic system, influenced by several genetic and environmental factors which may alter the course of the infection and its reproducibility. The inoculation method proved to be critical for the assessment of virulence alterations, with the drop inoculation on detached leaves resulting in the most stable assay, and the spray on whole seedlings resulting in the most extreme phenotypes, with ∆modcl2 and ∆modcl1/2 showing a complete loss of pathogenicity. These results both confirm and expand the key role of MoDCL2 first hypothesized by Raman and colleagues in a variety of fungal biological processes, including: sRNA biogenesis, fungal development and, as shown here, fungal virulence [27,50]. The reduction of virulence observed among RNAi mutants throughout the different setups further substantiate the importance of a fully functioning fungal silencing machinery and sRNAs generated for a successful infection of Bd tissues. GOE analysis of Mo DEGs highlighted an enrichment of terms mainly related to fungal growth and development in all infected tissues. Especially in the root GO, terms related to interactions with other organisms, particularly via protein secretion, were found overrepresented. By assessing expression of protein effectors, we confirmed upregulation of known avirulence genes, including AVR Pita1 (coding for a metalloprotease) and PWL2 (a glycine-rich protein recognized by weeping lovegrass and finger millet R proteins), both found highly upregulated in the foliar biotrophic phase (2 DPI), and Avr-Pik/km/kp, upregulated at 2 and 4 DPI.
Genes involved in appressorium and penetration peg functionality were upregulated in all three setups, including WISH (Water wettability, Infection, Surface sensing and Hyper-conidiation) G-protein coupled receptor protein, whose KO renders the fungus unable to develop appressoria and establishes the infection on intact rice leaves [54] and superoxide-generating NADPH oxidases NOX1 (MGG_00750) and NoxD (MGG_09956). NOX1 is involved in cell wall biogenesis, affecting both chitin and melanin biosynthesis and deposition. KO of this gene resulted in a complete loss of appressorium-mediated cuticle penetration and failure of in planta proliferation even when inoculated onto wounded rice seedlings [55]. Similar defects were detected in rice leaves and roots inoculated with ∆noxD [56], confirming the key role played by NOX proteins not only in growth and sexual reproduction, but also in fungal virulence. PLC3, significantly upregulated at both timepoints in leaves, is also required for normal conidiation and appressorium function, but contrary to NOXs normal infection was observed when ∆moplc3 was inoculated by infiltration into wounded rice leaves, excluding a function of PLC3 in fungal growth once inside host cells [57]. As expected, we detected in the later timepoint of the leaf infection (4 DPI) upregulation of genes known to be required for infection maintenance and expansion. One of these genes, S-(hydroxymethyl) glutathione dehydrogenase (MGG_06011), was exclusively upregulated in the 4 DPI dataset and it is known to be involved in the growth of infectious hyphae on barley leaves [58].
Overall, our analysis shows that Mo undergoes an extensive reprogramming during the establishment and maintenance of the infection and highlights the commonalities and differences in expression patterns depending on the inoculated tissue and the progression of the colonization. Of note and consistent with the detection of necrotrophic lesions on Bd roots ( Figure S4), fungal reprogramming during the root infection is closer to the one reported in the foliar necrotrophic phase (with 701 DEGs shared between the setups) compared to the biotrophic phase (only 392 DEGs shared between the setups).
mRNA sequencing of Mo-infected and non-infected Bd leaves and root samples allowed for a systemic analysis and comparison of expression pattern alterations of Bd genes in response to the fungal pathogen. As predicted, the highest amount of DEGs was observed at 4 DPI in leaves, when the infection is spreading outside of the inoculation point and the fungus has switched to a more aggressive necrotrophic lifestyle [5,12,59]. The relatively low amount of DEGs in the other two setups is consistent with the early infection stages and the limited amount of fungal biomass both in the root and at 2 DPI in leaves. Interestingly, in all three setups the percentage of upregulated DEGs was higher compared to downregulated DEGS, with 96% upregulated in leaves at 2 DPI, 60.7% at 4 DPI and 69.5% in the roots, indicating an overall strong induction of gene expression. Consistent with the highest numbers of DEGs detected at 4 DPI, most GO terms reported in Table S5 belong to this dataset, with terms related to metabolic and biosynthetic processes being the most prevalent. Interestingly, all datasets had enriched terms related to oxidoreductase activities (GO:001649) confirmed also by the upregulation of BdiBd21-3.4G0171000, coding for a multicopper oxidase, in all setups, and BdiBd21-3.1G0233800 and BdiBd21-3.1G0233900, coding for peroxidases, in the leaf 4 DPI datasets. Peroxidases are commonly associated with plant responses to stress and specifically fungal infection, as they are involved in a variety of processes including the synthesis of phenols such as tannins and melanins, reactive oxygen species (ROS) removal, lignin biosynthesis and induction of defense responses by stimulating intracellular Ca2+ signaling [60]. Another gene upregulated in all three setups is secologanin synthase-like (BdiBd21-3.2G0563800), encoding an enzyme involved in the biosynthesis of monoterpenoid indole alkaloids (MIAs), also reported upregulated in Bd following F. pseudograminearum infection [61]. Expression of other defense-related genes was induced in one or more of our datasets, including known pathogenesis-related proteins, receptor kinases, transcription factors and ABC transporters. Once again, 2 DPI and root samples were found overlapping in only a few of these genes. Examples of genes consistently found upregulated in all three setups are BdiBd21-3.3G0144800, encoding a protein kinase xa21 which confers resistance to Xanthomonas oryzae pv. oryzae race 6 [62], an ABC transporter A family member 2-like (BdiBd21-3.3G0465100) and the pleiotropic drug resistance protein 3-like (BdiBd21-3.2G0550500).
Interestingly, all PR genes shown in Table 5 were upregulated in leaves at 4 DPI, while leaf 2 DPI and root PR protein expression patterns did not overlap, with only two (PR-like and PR1-like) upregulated at 2 DPI, and three (two PR1-like and PR10) in the root datasets. Similarly, the strongest and most widespread upregulation in the other pathogen sensing/defense-related genes was observed at 4 DPI, with an upregulation of transcription factors belonging to MYB, WRKY and NAC families consistent with the upregulation observed in Bd after F. pseudograminearum infection [61], and the knowledge that these families regulate a variety of plant responses, including those to biotic stresses [63][64][65]. In line with the considerable overlap in Mo gene expression patterns between leaf 4 DPI and root datasets, the majority of genes found upregulated in Bd roots (80.9%) are also detected upregulated in leaves at 4 DPI, compared to 18% shared with the 2 DPI sample. Altogether, these results highlight differences in the upregulation of specific protein family members depending on the infected tissue and fungal lifestyle, while confirming the relevance of these families in response to Mo infections.
Candidate Bd miRNA-generating loci were identified with Shortstack from both control and infected filtered sRNA datasets. Given that the program-assigned cluster IDs are specific to the dataset analyzed and are not comparable between different files, we compared clusters based on their genomic coordinates, and further analyzed loci from the infected sample that had higher RPM (reads per million) than in control, or that were found exclusively in the presence of Mo. This resulted in the identification of 12 miRNA stem loop precursor sequences and 17 mature miRNAs across setups.
Specifically, we identified five mature miRNAs belonging to the MIR156 family, shown to be induced by F. oxysporum in Persicaria minor [66] and known to be induced by environmental stress, resulting in the cleavage of a SPL (SQUAMOSA promoter binding proteinlike) transcription factor and overall modulation of anthocyanin biosynthesis [34,67]. Bd-MIR529, also predicted to target this gene, was found upregulated at 4 DPI in our datasets. Another miRNA involved in the plant response to abiotic stress is miR159b, found upregulated during Mo leaf infection both at 2 and 4 DPI, and known to target two MYB transcription factors in cucumber (CsGAMYB1 and CsMYB29-like), involved in ABA signaling [68]. Finally, four additional miRNA families were detected as induced by Mo infection (MIR9484, miR9481b, MIR531 and miR7723a).
To establish the origin of the sRNA reads detected in the different leaf and root setups of the Mo--Bd interaction, sRNAseq datasets from infected samples were aligned to both the Bd 21-3 and the Mo 70-15 genome. Only reads aligning without mismatches to Mo and with at least two mismatches to Bd were assigned to the fungus and vice-versa, only reads aligning without mismatches to Bd and with at least two mismatches to Mo were assigned to the plant. As expected from the low amount of Mo in infected samples from leaves at 2 DPI, most of the reads were assigned to Bd, whereas higher levels of Mo reads were detected at 4 DPI (leaf) consistent with proliferating infection. All assigned reads were then filtered based on their read counts to select only reads either induced or upregulated in the datasets of infected tissues compared to uninfected tissues and axenic mycelia. We noted that most of the reads (>50%) found in infected samples are specific and are not detected in healthy tissues and axenic culture, showing that sRNA production both in the plant and the fungus is strongly responsive to infection. It follows that sRNA datasets from healthy plants and axenic fungal culture do not represent the full diversity of infection-related sRNA communities. As an additional step, we selected for sRNAs that were not aligning to the coding sequences of the organism of origin. The reasoning behind this filtering step is that we avoided accidental mRNA degradation to be kept as candidate ck-sRNAs, and more important, removed the sRNA sequences more likely to play an endogenous role. Given that the size distribution of upregulated/induced sRNA reads did not show variation in peaks compared to the total sRNA reads, we decided to select 20-22 nt sRNAs (canonical length range for PTGS) for further analysis. Target prediction was carried out with psRNATarget, a web-based prediction software specifically designed for plant sRNA investigations, which allowed for the identification of complementary mRNA sequences in the interacting organism. In PTGS, sRNAs are loaded onto AGO proteins, which guide them towards a complementary mRNA sequence that will then be degraded or sequestered, resulting in reduced levels of the encoded protein. Knowing the expression levels of the predicted targets from the same biological samples used for the sRNA sequencing, we proceeded with further selection of candidate ck-sRNAs based on the significant downregulation of their mRNA targets. Most of the predicted Mo sRNA effectors in the 2 DPI leaf (biotrophic phase) and root samples did not pass this filtering step, as their predicted targets were either upregulated or had comparable expression levels in the corresponding control datasets. There are a few possible explanations as to why the potential targets would not be significantly downregulated in our mRNAseq datasets, including: (i.) the sRNA has not yet been transported throughout the tissue, so the downregulation occurs only at the penetration site, where the fungus is physically interacting with the plant cells, and that is masked by the upregulation in distal parts of the tissue; (ii.) the target mRNA is not cleaved, but its translation is inhibited by the RISC complex acting as a physical barrier, in which case the measurable effect would not be at the mRNA level but only at the protein level; and (iii.) the target is indeed cleaved, but concurrently with the downregulating effect of the sRNA, there is a stronger endogenous upregulation of the gene, leading to either similar levels of mRNA as the control, or even higher.
Comparison of lists with predicted fungal mRNA targets of Bd-sRNAs between the different setups highlights substantial target conservation between the leaf biotrophic and necrotrophic phases, with 140 predicted shared Mo targets between the two, and 42 Mo targets conserved among all three setups. Subjecting the shortlisted predicted fungal target genes to a PHI-base survey for mutations in Mo with lethal or detrimental outcome, we found a clear indication for a reduction in pathogenicity or loss of virulence in respective KO mutants, including MoATG17 (MGG_07667), an autophagy-related protein whose KO impaired appressorium formation and function, resulting in a complete lack of disease symptoms on rice leaves [45]. Similarly, MoATG1 (MGG_06393) was also predicted to be targeted by Bd sRNAs and downregulated at 2 DPI, with Mgatg1 mutant reported by Liu et al. to be unable to infect rice and barley leaves [69]. Additionally, we predicted sRNAs targeting the avirulence effector AvrPiz-t (MGG_09055). AvrPiz-t suppresses rice PTI signaling pathway by targeting the E3 ubiquitin ligase APIP6 and suppressing its ligase activity, resulting in reduced flg22-induced ROS generation and overall enhanced susceptibility in vivo [70].
Interestingly, we detected Calpain-9 to be targeted in all three datasets and significantly downregulated in the foliar ones, both at 2 and 4 DPI. While downregulation in the root setup did not reach a significant padj value, targeting of this transcript would match the finding in the cotton-V. dahliae pathosystem, where the host was found to be expressing and exporting miRNAs to the infecting fungus to inhibit its virulence via the targeting of Clp-1 [39].
Additional predicted targets included integral membrane proteins (MGG_04378T0, MGG_04935T0), ergosterol biosynthesis enzymes (ERG6, MGG_10568T0; ERG26, MGG_ 04938T0) and fungal cell wall related genes, such as GPI-anchored cell wall beta-1,3-endoglucanase EglC (MGG_10400T0, targeted and significantly downregulated in all three datasets), cell wall protein MGG_09460T0 (targeted and significantly downregulated in both foliar samples), chitin deacetylase (MGG_08774T0) and chitinase 1 (MGG_08054T0). Targeting fungal membrane components and ergosterol homeostasis has already been proven to be a successful strategy in crop protection against fungal pathogens, both via the application of DMI (sterol demethylation inhibitors) fungicides [71] and with the RNAi-based HIGS (Host-Induced Gene Silencing) and SIGS (Spray Induced Gene Silencing) approaches, where artificial sRNAs are introduced in the plant either via transformation or external application is then transferred to the fungal cells during infection [72][73][74]. It is interesting to observe how the plant appears to have naturally evolved to produce sRNAs potentially able to target fungal essential and pathogenicity related genes. However, Mo is still able to progress its infection both in Bd leaves and roots, and the most likely factors behind this fungal success are the countermeasures it employs in this crosstalk, from the extensive reprogramming of gene expression to the release of effectors and sRNAs.
To substantiate the hypothesis that fungal sRNAs function as effectors to aid the establishment and maintenance of infection, we investigated the role of downregulated Bd targets. Targeting conserved sequences, such as ribosome-and photosynthesis-related ones, and hampering gene expression and biosynthetic processes would prove to be a more effective strategy than specific defense/immunity genes, which are more prone to mutate in the arms race between plants and pathogens.
Altogether, the prediction of Mo ck-sRNAs and corresponding Bd targets involved in the plant response to biotic and abiotic stress paves the way for further validation of predicted sRNA/mRNA interactions, including: (i.) proof of target cleavage via RACE or degradome sequencing; (ii.) verification of sRNA-target interaction in transient expression systems such as leaves of Nicotiana benthamiana; (iii.) mutational KO strategies of predicted target genes and/or precursor loci of predicted ckRNAs; and (iv.) detection of direct association of sRNAs or their target mRNA with the respective fungal or plant AGO1-like protein by immunopurification techniques [79,80].

AGO and DCL Protein Analysis and 3D Structure Modeling
Known ARGONAUTE and DICER-like protein sequences were downloaded from the NCBI database and analyzed following the workflow utilized for B. distachyon AGOs and DCLs [35]. The phylogenetic analysis and tree rendering were done by the Phylogeny.fr web server [81]. Multiple sequence alignment (MSA) of AGO PIWI domains was done using Clustal Omega [82,83] and visualized with the Mview multiple alignment viewer [84]. Protein domains were identified by using Simple Modular Architecture Research Tool (SMART) including PFAM domains in the search [85,86] and visualized with the Illustrator for Biological Sequences (IBS) online illustrator [87]. Prediction of protein location was done using the plant subcellular localization integrative predictor (PSI) [88], and prediction of the interactome was done using the STRING database [89], excluding text mining as indication of putative interaction. Finally, AGOs and DCLs amino acid (aa) sequences were utilized for predicting the proteins' 3D structure utilizing SWISS-MODEL [90]. Predictions were selected for further validations based on the GMQE and QMEAN Z-score values [91]. PROCHECK [92,93] and WHATCHECK [94] were used to check the stereochemical quality of the selected structures and calculate the Ramachandran Z-score [95]. Open-Source Py-MOL (The Py-MOL Molecular Graphics System Version 2.4.0a0) was used for visualization of the predicted structures [96].

Mo Mutants Cultivation and Inoculation
The Magnaporthe oryzae wild type Mo 70-15 and mutants ∆moago1, ∆moago2, ∆moago3, ∆modcl1, ∆modcl2 and ∆modcl1/2 obtained from N. Donofrio, Newark, NJ, USA, were grown as described [12]. Conidia germination and appressoria development was assessed by incubating conidial suspension (2 × 103 spores mL -1 ) in distilled water on poly-Llysine coated glass slides (Sigma-Aldrich, St. Louis, MO, USA) in a damp chamber at room temperature for 24 h and examined via inverted microscopy. Fungal stock was prepared on oatmeal agar (OMA) in a regime of 26 • C/24 • C (day/night cycle) and a light intensity of 70 µmol m −2 s −1 photon flux density. Seeds of Brachypodium distachyon cv. 'Bd21-3 [97] were germinated in soil (Fruhstorfer Erde Typ T) and cultivated in a growth chamber at 22 • C/18 • C (day/night cycle) with 60% relative humidity and a photoperiod of 240 µmol m −2 s −1 photon flux density. Three methods were utilized to assess disease progression and phenotype of the mutants on leaves: (i.) spray infection of whole seedlings on three-week-old Bd plants with Mo conidia suspension of 120 × 103 spores mL −1 in 0.002% Tween20 and assayed on the second youngest leaf; (ii.) spray infection of second youngest detached leaves of three-week-old Bd plants with 250 µL Magnaporthe oryzae conidia suspension 50 × 103 spores mL −1 in 0.002% Tween20; (iii.) drop inoculation on second youngest detached leaves of three-week-old Bd with 10 µL of conidia solution (50,000 conidia/mL in 0.002% Tween water) on 1% agar plates. Control leaves/seedlings were mock-inoculated with 0.002% Tween water in all setups. Leaves were kept at 16 h light (160 µmol m −2 s −1 )/8 h dark cycle at 22 • C/18 • C. Score disease progression and analysis of the necrotic spots was assayed 5 DPI via ImageJ software [98]. For root inoculation, sterilized seeds (3% NaClO for 15 min, followed by three times 5 min washes in sterile water) of Bd21-3 were vernalized in the dark at 4 • C for two days on half strength MS medium [99] and then moved to a 16 h light (160 µmol m −2 s −1 )/8 h dark cycle at 22 • C/18 • C. Roots of one-week-old seedlings were dip-inoculated in 1 mL of conidia solution (125,000 conidia/mL in 0.002% Tween water) for 3 h and transplanted in a (2:1) mixture of vermiculite (Deutsche Vermiculite GmbH, Sprockhövel, Germany) and Oil-Dri (Damolin, Mettmann, Germany). Control roots were mock-inoculated with 1 mL of 0.002% Tween water solution. Quantification of Mo DNA presence in roots was performed at 5 DPI by quantitative PCR based on the fungal actin (MoActin). The experiment was repeated two times, each time with n = 10 roots per experimental group.

Sample Preparation from Mo-Bd Interaction Sequencing
Mo was grown on oatmeal agar (OMA) for two weeks at 26 • C with 16 h light/8 h dark cycles both for sampling of mycelium and conidia production. Samples from axenic cultures were collected by scraping a mixture of mycelia and spores from three plates, followed by immediate freezing in liquid nitrogen. For root inoculation, sterilized seeds (3% NaClO for 15 min, followed by three times 5 min washes in sterile water) of Bd21-3 were vernalized in the dark at 4 • C for two days on half strength MS medium and then moved to a 16 h light (160 µmol m −2 s −1 )/8 h dark cycle at 22 • C/18 • C. Roots of oneweek-old seedlings were dip-inoculated in 1 mL of conidia solution (250,000 conidia/mL in 0.002% Tween water) for 3 h, transplanted in a (2:1) mixture of vermiculite (Deutsche Vermiculite GmbH, Sprockhövel, Germany) and Oil-Dri (Damolin, Mettmann, Germany) and grown for an additional 4 days before harvesting. Control roots were mock-inoculated with 1 mL of Tween water solution. For leaf inoculation, third leaves of three-week-old Bd were detached and drop-inoculated with 10 µL of conidia solution (50,000 conidia/mL in 0.002% Tween water) on 1% agar plates. Control leaves were mock-inoculated with Tween water. Leaves were kept at 6 h light (160 µmol m −2 s −1 )/8 h dark cycle at 22 • C/18 • C and collected for sequencing at 2 days post inoculation (DPI) and 4 DPI.

RNA Extraction, Library Preparation and Sequencing
Three roots or two leaves, respectively, were pooled per sample for RNA extraction and for each condition three pooled biological samples were prepared. Frozen tissue stored at −80 • C was ground in liquid nitrogen using mortar and pestle. Total

Transcriptome Analysis
Paired end sequenced cDNA reads of Illumina TruSeq ® Stranded mRNA libraries were analyzed through the quality check in FastQC and alignment in the junction mapper HISAT2 [100]. Magnaporthe oryzae MG8 release 38 [101] and Brachypodium distachyon Bd21-3 v1.1 (DOE-JGI, http://phytozome.jgi.doe.gov/) assemblies were utilized throughout this study as references. Htseq-count [102] and DESeq2 [103] were then used for read counting and differential gene expression calling (DGE) between the infected and control sample genes and to generate volcano plots. Heatmaps for selected DEGs were obtained with the pheatmap package for R [104]. Gene Ontology Enrichment (GOE) analysis on DEGs was done with AgriGO v2 [105]. Gene descriptions were integrated from the organism genome assembly, ENSEMBL Biomart, Phytozome and Blast2GO [106].

sRNA Analysis, Prediction of Endo-and Cross-Kingdom sRNA
The single end sequenced cDNA reads of Illumina TruSeq ® Small RNA libraries were analyzed starting with quality check with FastQC [107] and trimming of adapter artifacts with cutadapt [108]. The alignment of the reads to reference genomes and transcriptomes of Bd and Mo was done using the short read aligner Bowtie [109]. Reads with a 100% alignment to the genome of the organism of origin were selected, alongside the reads with at least two mismatches in the alignment to the target organism genome. Venn diagrams for sRNA and target overlaps were obtained with the VennDiagram package for R [110].
To identify interaction-related Bd sRNAs with endogenous function, both infected and control datasets were analyzed with ShortStack [111] to identify potential miRNA generating loci. Genomic coordinates and corresponding reads per million (RPM) of the identified clusters were compared between infected and control datasets to select clusters exclusively present or increased during infection. Both potential precursors and mature miRNAs deriving from these clusters were compared to known miRNA sequences, obtained from miRBase [112]. The structure of miRNA generating clusters was visualized with strucVis (version 0.4, Michael J. Axtell).
Bioinformatics analysis of candidate ck-sRNAs was done as described in Zanini et al. [43]. Only sRNA reads of 20-22 nt originating from non-coding regions and with a higher count in the organism of origin control datasets compared to the infected ones were analyzed further for ck-sRNA effector identification by the target prediction software psRNATarget used with customized settings [113].
Expression levels obtained for each gene were used as confirmation of downregulation of predicted targets from the psRNATarget software. PHI-base, a collection of experimentally verified pathogenicity/virulence genes from fungal and microbial pathogens [114], was used to gather information regarding phenotype and virulence of fungal mutants carrying a mutation in the identified Mo gene targets.

Conclusions
In the present work, we analyzed and characterized the interaction between Brachypodium distachyon and Magnaporthe oryzae at different fungal lifestyles and infection sites, both from a transcriptomic and sRNA expression profiles' point of view. The pathosystem has been studied as a model for the effect of blast disease on staple crops leaves (e.g., rice, wheat and barley), owing to Bd's short maturation phase, smaller genome and space-saving production [11,12,14]. In addition to foliar infections, we also established and characterized the interaction and responses to Bd root colonization by Mo. Additional to the confirmation of the extensive reprogramming in both organisms throughout the interaction, our results support the possibility that major staple crops co-evolved mechanisms of RNA-based communication with their microbial pathogens. Based on concomitant deep sequencing of mRNA and sRNA fractions, our work provides the first indication of both plant and fungal sRNAs involvement in the communication between Magnaporthe oryzae with the model grass Brachypodium distachyon, further supporting the theory of ckRNAi participation in plant-pathogen interactions. Interestingly, sRNAs induced during infection setups show only partial overlap both among the different tissues (leaves, roots) and the different infection phases (leaf: biotrophic, necrotrophic), raising the possibility that ckRNAi in a given host-pathogen interaction exhibits tissue-and lifestyle-specificity.
Supplementary Materials: The following are available online at https://www.mdpi.com/1422-006 7/22/2/650/s1, Figure S1: Analysis of MoAGO and MoDCL protein sequences, Figure S2: Multiple sequence alignment of the PIWI domain of MoAGO proteins, Figure S3: Development of appressoria from conidia of Mo wt and RNAi mutants, Figure S4: Phenotypic analysis of Mo wt and RNA interference mutants, Figure S5: Heatmap for Mo and Bd DEG calling with DESeq2, Figure S6: Results of gene ontology enrichment (GOE) analysis for significantly DE Bd genes in the 4 DPI leaf setup, Figure S7: Results of gene ontology enrichment (GOE) analysis for significantly DE Mo genes in the root setup, Figure S8: Visual representation of the identified upregulated Bd clusters (miRNA precursors) structures,