A novel plant-fungal association reveals fundamental sRNA and gene expression reprogramming at the onset of symbiosis


 Beneficial associations between plants and microbes are widespread in nature and have been studied extensively in the microbial-dominant environment of the rhizosphere. Such associations are highly advantageous for the organisms involved, benefiting soil microbes by providing them access to plant metabolites, while plant growth and development are enhanced through the promotion of nutrient uptake and/or protection against (a)biotic stresses. While the establishment and maintenance of mutualistic associations have been shown to require genetic and epigenetic reprogramming, as well as an exchange of effector molecules between microbes and plants, whether short RNAs are able to effect such changes is currently unknown. Here, we established an interaction between the model grass species Brachypodium distachyon (Bd, Pooideae) and the beneficial fungal root endophyte Serendipita indica (Si, syn. Piriformospora indica, Sebacinales) to elucidate RNA interference-based regulatory changes in gene expression and small (s)RNA profiles that occurred during establishment of a Sebacinalean symbiosis.
 Colonization of Bd roots with Si resulted in higher grain yield, confirming the mutualistic character of this interaction. Resequencing of the Si genome using the Oxford Nanopore technique, followed by de novo assembly yielded in 57 contigs and 9441 predicted genes, including putative members of several families involved in sRNA production. Transcriptome analysis at an early stage of the mutualistic interaction identified 2963 differentially expressed genes (DEG) in Si and 317 in Bd line 21-3. The fungal DEGs were largely associated with carbohydrate metabolism, cell wall degradation, and nutrient uptake, while plant DEGs indicated modulation of (a)biotic stress responses and defense pathways. Additionally, 10% of the upregulated fungal DEGs encode candidate protein effectors, including six DELD proteins typical for Sebacinales. Analysis of the global changes in the sRNA profiles of both associated organisms revealed several putative endogenous plant sRNAs expressed during colonization belonging to known micro (mi)RNA families involved in growth and developmental regulation. Among Bd- and Si-generated sRNAs with putative functions in the interacting organism, we identified transcripts for proteins involved in circadian clock and flowering regulation as well as immunity as potential targets of fungal sRNAs, reflecting the beneficial activity of Si.
 We detected beneficial effects of Si colonization on Bd growth and development, and established a novel plant-mutualist interaction model between these organisms. Together, the changes in gene expression and identification of interaction-induced sRNAs in both organisms support sRNA-based regulation of defense responses and plant development in Bd, as well as nutrient acquisition and cell growth in Si. Our data suggests that a Sebacinalean symbiosis involves reciprocal sRNA targeting of genes during the interaction.



Conclusions:
We detected beneficial effects of Si colonization on Bd growth and development, and established a novel plant-mutualist interaction model between these organisms. Together, the changes in gene expression and identification of interaction-induced sRNAs in both organisms support sRNA-based regulation of defense responses and plant development in Bd, as well as nutrient acquisition and cell growth in Si. Our data suggests that a Sebacinalean symbiosis involves reciprocal sRNA targeting of genes during the interaction.
Keywords: Brachypodium distachyon, Genome sequencing, Sebacinalean symbiosis, Serendipita indica, Small RNAs Background Mutualistic associations between beneficial microbes and plants are widespread and highly advantageous, especially in the microbial-dominant environment of the rhizosphere. This relationship benefits soil microbes by providing them access to plant metabolites; in return, they enhance plant growth and development by promoting nutrient uptake and/or protection against (a)biotic stresses [1,2]. The beneficial or parasitic outcome of a plant-microbe interaction is governed by the genotype and physiological status of the host, identity of the microbe, and environmental factors such as soil type and nutrient availability [3,4]. The establishment and maintenance of mutualistic associations (here called symbiosis) require genetic and epigenetic reprogramming and metabolome modulation by the exchange of effector molecules between the beneficial microbe and the plant [5,6]. Beneficial microbes have a significant impact on crop production, due to their effects on plant health and yield. However, considerable gaps in knowledge prior to their establishment in agricultural practice remain, including systemic identification of microbial abundance and diversity in various ecosystems, understanding the influence of climate, soil conditions, management practices, and, lastly, elucidating the intricacies of molecular mechanisms governing establishment of colonization and nutrient acquisition [7].
Crucial for regulation of gene expression, RNA interference (RNAi) is a well-known eukaryotic gene silencing mechanism [8], mediated by small RNAs (sRNAs) of 20-24 nucleotides (nt) in size and RNAi-associated proteins, primarily from the Argonaute (AGO), Dicerlike (DCL) and RNA-dependent RNA polymerase (RdRP) families [9]. DCLs generate sRNAs from longer RNA molecules, whereas AGOs bind sRNAs within an RNA-induced silencing complex (RISC). In the context of plant-microbe interactions, microbial protein effectors are known to promote pathogen colonization by suppressing host immune responses [10] and have been described in mutualistic associations as well [5], including the Sebacinalean symbiosis [11]. Recent findings suggest that sRNAs, through RNAi-based regulatory mechanisms, also can serve as effectors of pathogenic microbes [12], whereby the sRNA is secreted to suppress translation of a host mRNA via RNAi. Conversely, plants can secrete sRNAs that target virulence-associated mRNAs in the microbe [13]. This transfer of sRNAs and subsequent gene silencing in the target organism is called cross-kingdom RNAi [12].
Si is an endophytic fungus belonging to the order Sebacinales that colonizes the rhizodermis and cortex of a large spectrum of plants [15]. Si serves as an excellent model for beneficial microbes as it (i) primes plants for disease resistance against biotrophic [16] and necrotrophic [17] fungi, oomycetes [18], and viruses [19]; (ii) enhances the tolerance of plants against abiotic stress [20]; (iii) promotes growth and yield [21]; (iv) has a sequenced 25 Mb genome [22]; and (v) is genetically transformable and culturable in axenic conditions [23]. Si initially undergoes a biotrophic growth phase during Arabidopsis (Arabidopsis thaliana) and barley (Hordeum vulgare) colonization, with suppression of innate immune responses [24,25] and activation of induced systemic resistance [16]. Subsequently, Si colonization of barley enters a cell-death associated phase and switches to a saprophytic lifestyle [26,27].
Bd is a temperate grass species belonging to the Pooideae subfamily and a model for genetic studies of stress resistance and yield parameters of cereals [28]. Bd is self-pollinating, genetically transformable, easy to cultivate, and has a sequenced genome of 272 Mb [29][30][31]. It shares evolutionary proximity and broad synteny with complex crop genomes, such as wheat and rice [14], and is a host for major cereal pathogens [32]. Additionally, RNAi is operational in Bd, with proven alteration of micro RNA (miRNA) expression patterns in response to abiotic stresses [33,34]. In silico analyses revealed that the Bd genome, similar to other cereals, contains an expansion of DCL and AGO families [35].
Currently, the significance of cross-kingdom RNAi in mutualistic interactions is largely unknown. A recent in silico study predicted that the arbuscular mycorrhizal fungus Rhizophagus irregularis generates sRNAs, which have predicted targets in the host plant Medicago truncatula [36]. Moreover, tRNA-derived sRNA fragments from rhizobial bacteria were shown to regulate host nodulation-associated genes by utilizing the host's RNAi machinery [37]. To investigate the role of sRNAs in another agronomically relevant mutualistic interaction, we established a protocol for Si colonization of the model Pooideae Bd. Additionally, integrative highthroughput sequencing and transcriptome analysis assessed symbiosis-associated changes in the mRNA and sRNA expression patterns of both organisms. We discuss here possible sRNA-based regulation that might be critical for the establishment of the Sebacinalean symbiosis.

Brachypodium distachyon Bd21-3 forms a mutualistic interaction with Serendipita indica
To investigate whether Bd can develop a beneficial interaction with Si, we established an inoculation protocol using one-week-old seedlings of Bd line Bd21-3, with dip-inoculation in 5 × 10 5 chlamydospores ml −1 for 3 h. Comparison of grain production in fully mature, colonized vs. non-colonized plants grown in soil showed that Si increased the number of filled grains/plant by 49.9% (Fig. 1a, Additional file 1), and total grain weight/plant increased by 38.1% (Fig. 1b, Additional file 1). Consistent with the observation that Si-colonized Bd21-3 plants flower several days earlier than control plants, they exhibited a 32.2% increase in the number of spikelets at 2 months after inoculation (Fig. 1c, Additional file 1). Concordantly, growth and biomass analyses of Bd21-3 seedlings revealed a significant 8.6% increase in shoot length (Fig. 1d) upon Si colonization.
Further analysis of Bd21-3 seedlings indicated that Si colonization increased lateral root growth, as early as 4 days post inoculation (4 DPI, Additional file 2: Figure  S1a). By 25 DPI, roots showed a more extensively branched structure (Additional file 2: Figure S1b). Microscopy of Si-inoculated Bd21-3 roots confirmed root surface colonization and proliferation of fungal spores after staining with chitin-specific WGA-AF 488 at 4 DPI ( Fig. 2a-d) and further on at 7 DPI and 14 DPI (Additional file 2: Figure S2). Inter-and intracellular colonization of Bd21-3 cells in the root differentiation zone also was visible after WGA-AF 488 and propidium iodide staining (Fig. 2e-h). These results suggest that establishment of a mutualistic symbiosis correlates with observable phenotypic changes by 4 DPI; thus, this time point was used to further investigate the Bd21-3-Si system.

Resequencing of the Si genome
To improve Si assembly, the genome was resequenced using MinION (25,167 reads, 324 Mb) and MiSeq (18,225,814 reads, 5.46 Gb); together they yielded approx. 6.0 Gb of sequence information. De novo assembly of the Nanopore sequence reads generated 57 contigs, accounting for a total length of 24.7 Mb and a N50 of 1.3 Mb. The draft genome sequence features GC content of 50.8%, similar to the first genome version with 2,359 contigs and GC content of 50.7% [22]. Analyses using the eukaryotic gene prediction tool Genemark-ES 4.33 [38] revealed 9441 predicted genes (75% of the genome; 59,045 exons), 9498 intergenic regions (25% of the genome), and a gene density of 380.68 genes/Mbp (Additional file 2: Table S1). Annotation of the Si genes using a GenDB version designed to process eukaryotic genomes possessing multi-exon genes [39,40] revealed that 4756 have a predicted function. Comparison of predicted genes from the resequenced Si genome (Si-2020) vs. the 2011 assembly [22] indicated that the vast majority are shared (90.3%), with 915 genes unique to the Si-2020 genome (Additional file 2: Figure S3). There is a reduction in gene model numbers relative to the 2011 assembly, which can be attributed to improved gene prediction tools for eukaryotic organisms and a considerably reduced number of contigs. Additionally, Si shares 2585 genes with another member of the Sebacinales, Serendipita vermifera, while 156 genes are shared only with the ectomycorrhizal fungus Laccaria bicolor and 2729 genes are common in all three species (Additional file 2: Figure S4).

Establishment of the Si-Bd symbiosis is associated with extensive transcriptional reprogramming
To assess how the symbiotic interaction affects gene expression in both organisms, mRNA was sequenced and analyzed (Additional file 2: Table S2) from the roots of Si-colonized Bd21-3 seedlings (sample Bd-Si) and mocktreated plants (Bd-C) at 4 DPI, and from 4-week-old axenic Si cultures (Si-ax). Comparison of reads between Bd-Si and Si-ax identified 2963 differentially expressed fungal genes (DEGs, DESeq2: Wald test, Benjamini-Hochberg (BH) adjustment, padj < 0.05), which accounts for 31.4% of the 9441 predicted Si genes. Comparison of reads from Bd-Si and Bd-C revealed 317 plant DEGs (0.66% out of approximately 47,917 protein-coding transcripts disclosed in the JGI v1.1 annotation, padj < 0.05). The interaction-responsive DEGs in Si and Bd21-3, split into up-and downregulated groups are shown in Fig. 3.
All significant DEGs were submitted to gene ontology term analysis against the reference background for Bd and a customized Si-specific background. The resulting enriched terms relate to metabolic, mainly redox processes and catalytic activity functions (Additional file 2: Table S3). Of the 25 highly downregulated Si DEGs, many encode proteins associated with metabolic reprogramming networks involved in nutrient exchange and adaptation to nutrient availability (Table 1). By contrast, highly upregulated Si DEGs encode for proteins involved in fungal catalytic and hydrolytic processes. This suggests that by 4 DPI, Si has entered a saprophytic-like growth phase similar to that detected in barley roots [22]. To investigate whether any Si DEG encodes effector proteins, a computational pipeline [41] was used to mine 982 genes identified in the Si-2020 genome that encode signal peptide-containing proteins, resulting in 480 putative protein effector genes. In total, 174 (36%) of these were significantly upregulated during colonization of Bd21-3 (Additional file 2: Table S4), including six DELD family proteins [22]. In Bd21-3, many of the highly downregulated interactionresponsive DEGs encode transcription factors or proteins associated with stress responses or circadian clock regulation. Those showing high levels of upregulation include genes linked to immune responses and hormone The results are from three independent biological replicates that are represented by green, black, and red dots on the scatter graph, individual data values for 1a-c in Additional File 1. For a and b, 1-week-old seedlings were inoculated with 5 × 10 5 chlamydospores per ml and grown for approximately 3 months in F-E type LD 80 soil; for c, spikelets were counted on 2-month-old plants; for d, seedlings were grown for 3 weeks on a vermiculite:oil dri mixture (semi-sterile conditions). The significance threshold for p values, corrected for multiple testing (Benjamini-Hochberg) was set at 0.05 (*≤ 0.05) and the Si effect was calculated as ((Mean_Si-Mean_Control)/Mean_Control) × 10 2 )  Imaging was done with a LEICA S8 confocal microscope (e-h: maximum projection; z-stack). For a-d, 1-week-old seedlings were inoculated with 5 × 10 5 chlamydospores per ml and subsequently grown on a plastic mesh over 0.5X MS; for e-h, Si-inoculated seedlings were grown on vermiculite:oil dri mix before harvesting at 4 DPI  signaling networks ( Table 2 and Additional file 2: Table  S3). In order to further validate our sequencing data, we confirmed the expression of five Si and five Bd21-3 DEGs from Tables 1 and 2 by RT-qPCR. Generally, the qPCR results show a similar fold change for the selected genes between the colonized root and the respective controls, compared to the mRNA-seq results (Additional file 2: Figure S5). Together, these results show that both organisms utilize a complex enzymatic arsenal to establish and control the symbiosis.

Prediction of Si RNAi genes
Since RNAi-mediated gene silencing has been documented in most but not all fungi [42], we assessed whether the Si-2020 genome encodes RNAi-related proteins with conserved domain architecture and homology to RNAi components in the model filamentous fungus Neurospora crassa [43]. Genes encoding predicted DCLs (G4U2H0, G4TBW9) with typical domains (dsRNA-binding, RNase III and helicase, [44]), QDE2-like proteins with PIWI domains typical of AGOs (G4TEK0, G4TLO4, [45]), an AGO-like protein (G4T5G9), and RdRPs (G4TNU7, G4TQP0) were identified and were expressed in axenic culture and Bd21-3-associated Si samples (Additional file 2: Table  S5). Thus, the Si genome is predicted to contain genes encoding critical components of the RNAi machinery. Based on these new data, and the earlier discovery that AGO and DCL families are expanded in the Bd genome [35], we decided to sequence the sRNAs of both organisms, in order to assess the role of RNAi-based regulation and communication in symbiosis.
sRNA profiles undergo a substantial change at the onset of the Si-Bd symbiosis To evaluate how the mutualistic interaction affects the sRNA profiles in the colonized root and respective Si and Bd controls, reads from Bd-C, Bd-Si, and Si-ax sRNA data sets were subjected to consecutive filtering steps (Additional file 2: Figure S6). This greatly reduced the number of raw reads to be analyzed and allowed us to distinguish between sRNAs with potential targets in the interacting organism (putative ck-sRNAs) and sRNAs with potential functions in the same organism (putative endogenous sRNAs) (Additional file 2: Table S6). Analysis of sRNAs from the Bd-Si dataset revealed that the total number of putative ck-sRNAs exceeds that of endogenous sRNAs in both Si (786,732 vs. 261,478) and Bd21-3 (17 million vs. 1.6 million), but the converse is true for unique sRNAs (

Size distribution profiles of Si and Bd21-3 sRNAs
Size distribution of sRNA reads from Bd21-3 and Si during colonization, and the respective controls was then assessed. For putative endogenous Bd21-3 sRNAs, peaks at 21 and 24 nt were identified, with the 24 nt sRNAs exhibiting greater diversity than those of 21 nt (Fig. 4a, b). These sizes are consistent with the expected peaks of RNAi-associated sRNAs in plants [46]. For putative endogenous Si sRNAs, a bimodal size distribution pattern was observed in the total fractions, with the first peak at 26 nt and second at 29-30 nt (Fig. 4c, d). A smaller peak of 21 nt long molecules was observed in the Bd-Si but not Si-ax samples, indicating that colonization affects the relative size distribution of Si sRNAs. Since previously identified ck-sRNAs range from 20 to 24 nt [12,13], the size distribution of putative ck-sRNAs and corresponding reads in the control samples was assessed over a narrower window. Contrary to endogenous sRNAs, ck-sRNAs displayed no prominent peaks in the 20-24 nt window (Additional file 2: Figure S7).
Before sRNAs can guide RNAi-mediated gene silencing, they must be loaded onto AGO proteins and assembled into a RISC. Previously, Arabidopsis AGO proteins were shown to preferentially recruit sRNAs with specific 5′ termini [47]. Hence, we analyzed the 5′ terminal nt composition of Bd-C, Si-ax, and Bd-Si sRNAs. For unique putative endogenous and ck-sRNAs, the 5′ nt composition was relatively consistent except for the 24 nt Bd21-3 sRNAs, which exhibited a strong bias towards 5′ A (Additional file 2: Figure S8, Figure S9). The total sRNA  fractions exhibited somewhat greater variability in 5′ nt composition. Of the total Bd21-3 endogenous sRNAs, 24 nt molecules from colonized and non-colonized tissue showed a strong bias towards 5′ A, while 21 nt molecules were biased towards a terminal U (Additional file 2: Figure  S10), and 20 nt ck-sRNAs had a higher percentage of 5′ Cs (Additional file 2: Figure S11). Of the total endogenous Si sRNAs, those from colonized samples generally had a stronger bias towards 5′ A than sRNA reads from Si-ax, especially at 26 nt and 21 nt (Additional file 2: Figure S10). A slightly higher percentage of 5′ As also was detected in total putative ck-sRNAs of 21 nt (Additional file 2: Figure S11).

Differentially expressed Si and Bd21-3 sRNAs
Analysis of unique plant sRNAs in Bd-C vs. Bd-Si revealed that 63% of the putative endogenous sRNAs were exclusively present in Bd-C, 30% were exclusively in Bd-Si and 7% were present in both (Fig. 5). For the reads from the putative ck-sRNA pipeline, 76% of the reads were exclusively present in Bd-C, 13% were exclusive to Bd-Si, and 11% were found in both. Comparison between the unique fungal sRNAs in Si-ax and Bd-Si indicated that 98.1% of the putative endogenous sRNAs were exclusively present in axenic culture, 0.98% were exclusively found in Bd-Si, and 0.92% were in both. Similarly, from the putative ck-sRNA pipeline, 98.2% of  the sRNA reads were exclusive to the Si-ax sample, 0.3 % were exclusive to Bd-Si, and 1.5% were present in both. Considering only the Si sRNAs in the Bd-Si sample, 51.1% of the putative endogenous ones and 15% of the putative ck-sRNAs are exclusively present in the colonized sample. Among Bd21-3 sRNAs in the Bd-Si sample, there are 80.5% putative endogenous sRNAs and 54.3% ck-sRNAs exclusive for the colonized root. These data show that colonization induces many novel putative endogenous and ck-sRNAs in Bd21-3, and a smaller amount in Si, due to fungal quantity in the colonized roots.

Identification of Bd21-3 miRNAs during the Bd-Si interaction
Using the ShortStack analysis tool, we identified Bd21-3 loci that correspond to putative endogenous sRNAs expressed during Si colonization. The DicerCall function indicated loci whose predominant sRNAs are 20-24 nt.
Comparison of these sRNAs with miRBase identified 16 sRNA-generating loci that correlate to known miRNAs (Table 3). These miRNAs belong to highly conserved plant miRNA families that regulate growth and development [33,34]. We conducted the same analysis with Si sRNAs, but no predicted miRNA-like RNAs were identified in the colonized sample, possibly due to a lack of data about the fungal sRNA-generating loci.
To assess whether Bd21-3 generates ck-sRNAs that potentially target Si genes, we searched for predicted targets of 329 Bd21-3 sRNAs (21 nt long) induced in Si-    Table S7). In order to confirm the expression of some of these sRNAs, we conducted stem-loop PCR and sRNA-specific sequencing on 10 Si and 10 Bd21-3-originating sRNAs from our Bd-Si sample and Table S7. All Si and all Bd21-3 sRNAs were amplified in the stem-loop PCR. To verify the nature of the amplification products, a subset of four SisRNAs and three BdsRNAs were cloned and sent for Sanger sequencing, confirming the expected sRNA sequences (Additional file 2: Figure S12, original gel pictures in Additional file 3 and Additional File 4, sequencing results in Additional file 2: Table S8). Thus, predicted targets of putative ck-sRNAs within this system imply another layer of expression control within the mutualistic interaction.

Discussion
We established and studied the interaction between Brachypodium distachyon-a model Pooideae plant with shared synteny to major cereal crops-and Serendipita indica-a beneficial endophyte with an exceptionally large host range. This particular combination of traits adorning the Bd-Si interaction has great translational value towards filling in the gaps in knowledge about plant symbioses, especially their transcriptomic and sRNA expression profiles and the significance of RNAi. We show that Si colonizes Bd, resulting in shoot growth promotion, earlier flowering, and improved grain development. In comparison, earlier studies characterizing sRNA expression was calculated as log 2 (colonized/control) from normalized reads, NA stands for sRNAs expressed exclusively in the colonized sample; transcript expression is indicated as the log 2 (colonized/control) FC the interaction between Si and barley have demonstrated that fungal hyphae establish an interface with the root cell plasma membrane at an early colonization stage, followed by expansion of an extracellular hyphal network, intercellular growth, and intracellular penetration of cortical and rhizodermal cells [26]. Around 3 to 5 DPI, Si starts the switch from a biotrophic to a saprophytic lifestyle [26,27]. Although this change involves intracellular hyphae extensively colonizing dead host cells and gradual digestion of cortical cell walls, the plant still benefits from the fungal presence. Consistent with these findings, our microscopic analyses confirmed proliferation of Si chlamydospores and both inter-and intracellular colonization of Bd21-3 cells in the root differentiation zone from 4 to 14 DPI. Detection of proliferating hyphae that were not wrapped in plasma membrane further suggests that Si is colonizing dead surface plant root cells at 4 DPI [25,26].

Transcriptional changes detected during the Bd-Si interaction
To investigate colonization of Bd21-3 by Si, we analyzed Si DEGs in colonized vs. axenic mycelium samples. Gene ontology analysis indicated enrichment in genes involved in various metabolic and catalytic processes. DEGs with the greatest changes in expression play roles in plant cell wall degradation, carbohydrate metabolism and nutrient acquisition. These changes in nutritional reprogramming are to be expected, considering the different nutritional content that Si has accessible in planta vs. axenic CM plates and a more detailed look into the roles of the changed genes unveils a typical switch of fungal lifestyle. Examples of upregulated Si genes involved in cell wall degradation include a probable Pectate lyase, Endo-1,4beta-xylanases, Cellulose 1,4-beta-cellobiosidase, and Rhamnogalacturonan acetylesterase. The genes encoding these hydrolytic enzymes, which have undergone expansion in the Si genome [49], are similar to those upregulated in Si during saprophytic growth on autoclaved barley roots at 3 DPI and 5 DPI [22]. PiAMT1, encoding a high affinity ammonium transporter also was upregulated (logFC = 3.35; padj < 0.0001). Other upregulated genes encode enzymes involved in carbohydrate metabolism, including probable glucosidases, glucanase, and galactosidase. These proteins may modulate glucose concentration, which then regulates expression of some cell-wall degrading enzymes [22,50]. Some of the 174 putative effector protein-encoding genes also are differentially expressed during Si colonization of barley and Arabidopsis [22,27]. Six of these proteins (Additional file 2: Table S4) contain the Si-specific DELD domain, which suggests that Si utilizes a common protein effector arsenal to colonize various hosts. Considering highly downregulated Si DEGs, several encoded proteins are associated with adaption to nutrient availability (Accumulation of dyads protein 2, ADY2; Succinate-fumarate transporter) and nutrient acquisition (Acyl-CoA dehydrogenase; Carnitine acetyltransferase, CRAT; Phenylalanine ammonia-lyase, PAL [51,52];). Their reduced expression suggests ample nutrient availability at 4 DPI. Given the similarities in the Si transcriptome during colonization of Bd and barley, we propose that these fungal-plant interactions follow a pattern, and that by 4 DPI, a network of plant-endophyte communication cues has initiated a tightly controlled transcriptional program, leading to a shift from biotrophic to saprophytic growth. Roots of Bd21-3 plants also displayed substantial transcriptional reprogramming following Si colonization. Gene ontology term analysis indicated enrichment in genes involved in catalytic and oxidoreductionassociated processes. Bd21-3 DEGs exhibiting the greatest changes in expression between colonized and noncolonized plants are related to stress-response, defense, and plant development. Of the downregulated Bd21-3 genes, several encode proteins commonly associated with stress responses, including a peroxidase, a woundinduced protein, and a putative protein kinase. Additionally, members of the Heat-shock protein gene family [53] are commonly induced in Bd during abiotic stress and members of the Abscisic acid/water deficit stress (ABA/ WDS)-induced protein and the Rare cold inducible (RCI2) gene families enhance abiotic stress tolerance in various plant species [54,55]. Circadian clock and flowering regulation genes such as Pseudo-response regulator 7 (PRR7), Cold regulated protein 27 and Constans-like protein (CO8), also are downregulated during Si colonization. While members of the PRR and CO protein families work together to control flowering time [56,57], any influence on early flowering in Si-colonized Bd21-3 is unclear. Circadian clock-associated genes also regulate lateral root development in Arabidopsis [58]; whether Si-induced changes in their expression influence root growth is unknown. Other downregulated development-associated DEGs include Fantastic four meristem regulator (FAF), which regulates shoot and root development [59], and putative Sulfoquinovosyltransferase (SQD2), which modulates seed setting and tiller development in rice [60]. Finally, several downregulated DEGs encode transcription factors, including MYB-related, GRAS, and bZIP.
In comparison, many of the upregulated Bd21-3 DEGs are associated with immune responses. Examples include genes encoding leucine-rich repeat (LRR) protein, a WRKY transcription factor, and thaumatin family protein. Increased expression of the defense gene Pathogenesis-related protein 1 (PR1) was similarly and transiently reported in Si-colonized Arabidopsis roots [61]. Upregulation of Glutathione S-transferase (GST) is consistent with the increased antioxidant capacity of Si-colonized plants, which provides protection against attack by necrotrophic pathogens [21,62]. The upregulation of genes in other hormonal networks (PGP-like Phosphoglycoprotein auxin transporter) and redox processes (Multicopper oxidase) further suggests that Si colonization affects a range of signaling pathways.

Bd miRNAs detected in the Bd-Si sample
The role of miRNAs as regulators of gene expression in the Sebacinalean symbiosis is largely unexplored. One report showed that Si induces growth promotionassociated miRNAs in Oncidium orchid roots [63]. Analysis of putative endogenous Bd21-3 sRNAs expressed during Si colonization identified 16 miRNAs. Some of them have known targets in transcription factors associated with plant growth and development. For example, the bdi-MIR166 family targets mRNAs encoding Homeobox domain-leucine zipper transcription factors [64]. In Arabidopsis, repression of these transcription factors by the miR165/166 family modulates root growth, maintenance of the shoot apical meristem, and development of leaf polarity [65]. Plant-specific transcription factors encoded by Squamosa promoter-binding protein-like (SPL) genes are the presumed targets of bdi-MIR156 and bdi-MIR529 [66]. In Arabidopsis, miR156-mediated downregulation of SPLs modulates developmental timing, lateral root development, branching, and leaf morphology [65]. Members of the MYB superfamily of transcription factors, which regulate many aspects of development, are the predicted targets of bdi-MIR159 [34]. Interestingly, miRNAs belonging to the miR159 and miR166 families in cotton are known ck-sRNAs that target virulence genes in Verticillium dahliae [13].
In combination with earlier studies on Bd RNAi proteins [35] and Bd interaction with the pathogen Magnaporthe oryzae [82], the in silico analyses presented here suggest that Si and Bd contain functional RNAi components and that both organisms generate ck-sRNAs, which potentially modulate this mutualistic interaction. However, further studies are necessary to validate crosskingdom RNAi in a Sebacinalean symbiosis. Namely, degradome analysis is needed to confirm target degradation and evidence that Bd21-3 and Si AGOs associate with sRNAs expressed by the interacting organism is necessary for confirmation of cross-kingdom RNAi.

Conclusions
We report that Bd21-3 and Si form a mutualistic symbiosis with a promoting effect on plant yield and development, accompanied by changes in gene expression in both organisms, including putative protein Si effectors and RNAi-related genes. sRNA profiles of both organisms also changed, indicating that this model system will provide important insights into the multiple layers of regulation and interaction between beneficial fungi and cereal hosts. Within the broader scope of plant-mutualist interactions, we show that detection of putative RNAiinvolved sRNAs in an interaction highly benefits from simultaneous transcriptome analysis and indicate an involvement of sRNA-based regulation in defense responses, nutritional reprogramming, and colonization maintenance. Alongside other experimental approaches in plant-microbe interactions (eg. sRNA uptake studies [83]), developing a deeper understanding of the communication mechanisms that modulate mutualistic interactions is highly relevant for establishing robust growth promotion and protection strategies in crops.

Bd and Si cultivation and inoculation
The seeds of Brachypodium distachyon (Bd) line Bd21-3 (gift from R. Sibout, INRA Versailles, France) were surface sterilized (3% active chlorine, sodium hypochlorite solution) for 15 min, washed three times, and placed on half-strength MS [84]  Samples for RNA-seq, RT-qPCR, stem-loop PCR, and microscopy were also grown under these conditions. To assess growth promotion in Si inoculated Bd relative to the control, we used the pairwise t test or the Mann-Whitney-Wilcoxon test on each of the three repetitions of experiments, after checking for normality and homogenous variances. Benjamini-Hochberg correction for multiple testing was used to correct the p values and the significance asterisks were assigned to the average p-value as follows: * for p ≤ 0.05, ** for p ≤ 0.001, and *** for p ≤ 0.0001.

Microscopy
Following Si inoculation, one-week-old seedlings were grown on plastic mesh (~90 μm) placed over halfstrength MS medium or on vermiculite/oil dri prior to assessing root colonization. Si was visualized with the chitin-specific dye WGA-AF 488 (wheat germ agglutinin; Molecular Probes, Karlsruhe, Germany), as described in Deshmukh et al. (2006) [26], with boiling in KOH (10%) for 30 s, prior to incubation in phosphate-buffered saline (PBS, pH 7.4). Root cells were visualized by incubating with propidium iodide (10 μg ml −1 ) for 10 min and washing with sterile water. Confocal laser scanning microscopy was done (TCS SP8 microscope, Leica, Bensheim, Germany) and the Leica LAS X software was utilized for visualization and maximum (z-stack) projections.
Resequencing, assembly, and annotation of the Si genome The MasterPure Yeast DNA Purification Kit (Epicentre, Illumina) was used to extract genomic DNA from 4-week-old axenic Si cultures. The Si genome was resequenced, assembled [86], and annotated as described [87], whereby a MinION sequencing library was prepared using the Nanopore Rapid DNA Sequencing kit. Sequencing was performed on an Oxford Nanopore MinION Mk1b sequencer using a R9.5 flow cell. Additionally, sequencing of an Illumina Nextra XT library was performed on the MiSeq platform (Illumina; 2 × 300 bp paired-end sequencing, v3 chemistry). Adapters and low-quality reads were removed by an in-house software pipeline prior to polishing [88]. MinKNOW (v1.13.1, Oxford Nanopore Technologies) was used to control the run with the 48 h sequencing run protocol, and base calling was performed offline using albacore (v2.3.1, Oxford Nanopore Technologies). The assembly was performed using Canu v1.6 ( [89], default settings). The resulting contigs were polished with Illumina short read data using Pilon [90] for eight iterative cycles. BWA-MEM [91] was used for read mapping in the first four iterations and Bowtie2 v2.3.2 [92] in the second set. Gene prediction was performed with GeneMark-ES 4.3.3. ( [38], default settings). Predicted genes were functionally annotated using a modified version of the genome annotation platform GenDB 2.0 [39] for eukaryotic genomes [40]. RNAi-associated proteins were predicted by searching the proteome [22] for typical domain structure and highest homology to Neurospora crassa RNAi proteins (NC12 genome assembly [93] [103]. An intron length of 20-2000 nt was allowed for Si [104] and 20-10,000 nt for Bd21-3 [105]. The reads were counted using HTSeq-count [106], differential gene expression was performed with DESeq2 [107], and gene enrichment analysis with AgriGO v.2 [108], with reference Bd21 setting for Bd21-3 (Bd 21 synonyms) and a customized background for Si. Volcano plots were generated using plotly [109] and ggplot2 [110] R [111] libraries. Gene descriptions were obtained from the organism annotations or Blast2GO [112].
sRNA analysis and prediction of putative endogenous and ck-sRNAs Raw reads from sRNA sequencing [113] were submitted to FastQC [101] and adapter trimming [114]. Bowtie [115] was used for alignment as detailed in Additional file 2: Fig. S6. The resequenced Si genome was used for alignments of fungal origin. tRNA/rRNA sequences were downloaded from RNAcentral ( [116], EMBL-EBI). Putative endogenous sRNA reads were submitted to Short-Stack [117]. For filtering putative ck-sRNAs, a previously established pipeline [118] was utilized. Reads were normalized to the total number of mapped reads for a single genome and reads per million (RPM) and log 2 (colonized/mock-treated) values calculated. Thus, a sRNA read was selected as a putative ck-sRNA if it was present exclusively or at a higher quantity (i.e., induced) in the colonized vs. control sample. Putative ck-sRNAs were submitted to psRNAtarget [119]. Since the separation of sRNA and mRNA fractions from each biological sample was facilitated, we checked for downregulation of mRNAs corresponding to predicted sRNA target genes within the DEGs. Transcriptomes used for these predictions were Bd 21-3 v1.1 (DOE-JGI, http://phytozome.jgi. doe.gov/ [102]) and Si DSM11827 ASM31354 v.1 [22].
For the identification of sRNAs in the interaction of Si with Bd21-3 stem-loop RT-PCR was employed [122]. cDNA was synthesized from DNase I-treated total RNA extracted from Si axenic culture or inoculated Bd21-3 roots. The folding of the hairpin primer was performed according to Kramer (2011) [123]. For each stem-loop reaction, six hairpin primers were multiplexed in a 20-μL reaction using the Revertaid RT enzyme according to the manufacturer's instructions (Thermo Scientific). For primer annealing, the reaction was incubated for 30 min at 16°C followed by an extension step at 42°C for 30 min. The reaction was stopped at 85°C for 5 min. cDNA was stored at − 80°C until further use. Endpoint PCR was performed using an universal stem-loop primer and specific sRNA primer (Additional file 2: Table S10) under the following conditions: initial denaturation at 95°C for 5 min followed by 35 cycles: 95°C for 30 s, primer annealing at 60°C for 30 s, and extension at 72°C for 30 s. PCR products were separated by gel electrophoresis on a 2% (w/v) agarose gel. To obtain sequence information of the amplified sRNAs of the stem-loop reaction, PCR products were purified and cloned into the pGEM®-T Easy Vector Systems (Promega, Madison, WI, USA) following the manufacturer's instructions. From each cloned sRNA, six colonies were further analyzed by Sanger sequencing using a M13 reverse primer.
Additional file 3. Uncropped gel picture annotated in Additional file 2: Figure S12a.
Additional file 4. Uncropped gel picture annotated in Additional file 2: Figure S12b.