Physiological stressors and invasive plant infections alter the small RNA transcriptome of the rice blast fungus, Magnaporthe oryzae

Background The rice blast fungus, Magnaporthe oryzae is a destructive pathogen of rice and other related crops, causing significant yield losses worldwide. Endogenous small RNAs (sRNAs), including small interfering RNAs (siRNAs) and microRNAs (miRNAs) are critical components of gene regulation in many eukaryotic organisms. Recently several new species of sRNAs have been identified in fungi. This fact along with the availability of genome sequence makes M. oryzae a compelling target for sRNA profiling. We have examined sRNA species and their biosynthetic genes in M. oryzae, and the degree to which these elements regulate fungal stress responses. To this end, we have characterized sRNAs under different physiological stress conditions, which had not yet been examined in this fungus. Results The resulting libraries are composed of more than 37 million total genome matched reads mapping to intergenic regions, coding sequences, retrotransposons, inverted, tandem, and other repeated regions of the genome with more than half of the small RNAs arising from intergenic regions. The 24 nucleotide (nt) size class of sRNAs was predominant. A comparison to transcriptional data of M. oryzae undergoing the same physiological stresses indicates that sRNAs play a role in transcriptional regulation for a small subset of genes. Support for this idea comes from generation and characterization of mutants putatively involved in sRNAs biogenesis; our results indicate that the deletion of Dicer-like genes and an RNA-Dependent RNA Polymerase gene increases the transcriptional regulation of this subset of genes, including one involved in virulence. Conclusions Various physiological stressors and in planta conditions alter the small RNA profile of the rice blast fungus. Characterization of sRNA biosynthetic mutants helps to clarify the role of sRNAs in transcriptional control.


Background
RNA silencing, the RNA-mediated suppression of a target gene in a sequence specific manner, has been described in many organisms as post-transcriptional gene silencing, RNA interference, and quelling, in plants, animals, and fungi, respectively [1,2]. These silencing pathways are involved in transposon silencing, viral defense, DNA elimination, heterochromatin formation, and post-transcriptional repression of cellular genes [3,4]. The critical steps of these pathways include production of small RNAs (sRNAs) of 19-27 nucleotide (nt) from structured or double-stranded RNA by RNase III-like endonucleases called Dicers, followed by loading of these sRNAs into Argonaute-containing complexes to form RNA-induced silencing complexes (RISC) that guide the cleavage of target transcripts [3,5].
Magnaporthe oryzae (hereafter M. oryzae), a filamentous fungal pathogen that causes rice blast disease, leads to severe yield losses in cultivated rice worldwide. The Magnaporthe-rice pathosystem is an excellent model to study the molecular basis of pathogenesis, in part because of the economic importance and molecular and genetic tractability of both host and pathogen, availability of genome sequences and transcriptome data, and the similarity of the rice blast fungus with other cereal diseases in terms of appressorium formation and intracellular tissue invasion [6][7][8]. The blast fungus infects rice plants at all stages and all tissuesleaves, stems, nodes, panicles and roots [8,9]. Leaf infection starts with reversible attachment of conidia to the host surface. The release of spore tip mucilage from conidia forms an irreversible attachment to the surface [7]. Conidia then germinate and form infection structures called appressoria at the end of the germ tubes. Appressoria generate enormous amounts of turgor pressure to penetrate the outer surface of the plant, which is followed by tissue invasion and colonization resulting in disease lesions [8].
The carefully controlled and timed life cycle of M. oryzae warrants examination into both transcriptional and post-transcriptional gene regulation. To date, endogenous small RNA species have been examined in several fungal species including the budding yeast, Saccharomyces cerevisiae, the model ascomycete fungus and close relative of M. oryzae, Neurospora crassa and M. oryzae itself ( [5]; reviewed in [10][11][12]). The M. oryzae genome encodes putative RNAi pathway components consisting of two Dicer-like proteins (DCL), three Argonautes (AGO), three RNA-dependent RNA polymerases (RdRP) and one exportin-5 [11,13,14]. In 2004, Kadotani et al. showed that Magnaporthe Dicer-like-2 (MoDcl2) is necessary for generating siRNAs in hairpin dsRNA-mediated RNA silencing, unlike their counterparts in Neurospora crassa in which both of the Dicers are redundantly involved in the silencing process [15]. While the Kadotani et al. study (2004) did not find a clear role for the Magnaporthe Dicer-like-1 (MoDcl1) in the generation of siRNAs, the function of this gene is yet to be discovered. Exportin-5, which plays a role in nuclear export of pre-miRNAs in animals, was identified and deleted in M. oryzae and shown to have reduced pathogenicity on roots, and to complement its yeast homolog, the nucleocytoplasmic transporter Msn5 [16][17][18].
The complexity of RNAi components in this fungus suggests the possible existence of multiple sRNA pathways. Though miRNAs have not been identified in M. oryzae, the published data have demonstrated numerous components of the M. oryzae sRNA repertoire, including small interfering (siRNA) -induced silencing of transgenes [11], the presence of viral derived siRNAs [19], methylguanosine-capped and polyadenylated small RNAs (CPA-sRNAs) [20] and sRNAs matching to repetitive elements, tRNA loci, rRNAs, protein coding genes, snRNAs and intergenic regions ( [12]; this study). Moreover, several new classes of small RNAs have been described recently in fungi, including QDE-2 -interacting siRNAs (qiRNAs), microRNA-like small RNAs (milRNAs) and Dicer-independent small interfering RNAs (disiRNAs) of N. crassa [21]. The existence of RNAi has been demonstrated in many fungal species including pathogenic species [22][23][24]. In Mucor circinelloides, production of two classes of small antisense RNAs requires a functional Dcl-2 gene [25].
In the M. oryzae life cycle, transition from the biotrophic to the necrotrophic stage is likely induced by environmental cues that the fungus encounters during the pathogenic process [26]. Conditions encountered by the fungus during its growth in rice leaves, particularly during disease symptom expression, include nitrogen starvation, and this is regarded as one of the inductive cues for disease symptom expression during rice plant infections [26,27]. Starvation is also considered a key factor for inducing the germ tube tip to differentiate into an appressorium [28]. These and similar studies prompted us to examine whether posttranscriptional changes occur in M. oryzae during exposure to different physiological stressors, as well as during several invasive stages of pathogenic growth. Recent studies on genome-wide transcriptional changes demonstrate that in vitro conditions such as nitrogen and carbon deprivation are comparable to genome-wide changes that occur during invasive plant growth [29]. The role of post-transcriptional control during these stresses, however, has not yet been studied. We examined small RNA profiles of M. oryzae during stress conditions matching those of Mathioni et al. [29], namely, nitrogen starvation, carbon starvation, minimal media, and paraquat-induced oxidative stress. We also examined small RNA fungal profiles in planta at 72 and 96 hours post inoculation (hpi), as well as mock inoculated plants. Small RNAs from certain conditions showed association with specific genomic loci; for example we observed an abundance of sRNAs associated with retrotransposons during nitrogen starvation. We also identified a small subset of genes whose transcriptional down-regulation during several of the stress conditions showed association with high numbers of sRNAs. To further examine whether this represented a true link between transcriptional and posttranscriptional control, we generated four small RNA biosynthetic fungal mutants, and performed quantitative real-time RT-PCR on this down-regulated subset of genes. These results revealed a previously unknown role for the Magnaporthe Dicer-Like1 (MoDcl1) gene.

Results
Small RNA profiling of M. oryzae stress libraries We analyzed small RNA (sRNA) size profiles of libraries constructed from M. oryzae mycelial tissue, and from susceptible rice leaves infected with M. oryzae at two time-points. For the former, M. oryzae was subjected to the in vitro stresses of nitrogen starvation (NS), carbon starvation (CS), minimal media (MM), and oxidativeinducing conditions via paraquat (PQ). Mycelia were also grown on a rich, complete media (CM) for comparative purposes. Growth conditions for these in vitro stresses mimicked those of Mathioni et al. [29], and hereafter will be referred to as the "mycelial" libraries. In order to compare the in vitro stress sRNA profiles to small RNAs generated during in planta growth, we also isolated fungal-induced rice lesions at 72 and 96 hpi. These will be hereafter referred to as the "in planta" libraries. Small RNA libraries were generated and quality checked before performing deep sequencing using Illumina's Genome Analyzer II (GA II ) system. This resulted in 58,511,280 and 28,941,539 sequences from the mycelial and in planta libraries, respectively ( Table 1). The sequences were aligned to the M. oryzae genome assembly version 6, downloaded from the Broad Institute (http://www.broadinstitute.org/). Total genomematched reads (mapping to one or more sites in the M. oryzae genome) and distinct genome-matched reads were calculated. There were 34,652,258 and 2,414,112 total genome-matched reads, and 927,501 and 54,394 distinct genome-matched reads for mycelial and in planta libraries, respectively, (Table 1). For the in planta libraries, M. oryzae sequences were parsed out from rice sequences by alignment to the fungal genome. An informatics workflow describing this process is presented in Additional file 1: Figure S1, along with an example of the output from the pipeline as applied to a mycelial library.
The size distribution of total genome-matched reads for mycelial libraries were bimodal with peaks at 23/24 nucleotides (nt) and 26/27 nt, while distinct genomematched reads showed a normal distribution curve with the peak at 23/24 nt ( Figure 1). The abundance of 24 nt reads represented 16.77% of the total population. In planta libraries also had a similar pattern of abundance, but the major peak was at 26 nt followed by 24 nt for the total genome-matched reads and 22/23 nt for the distinct genome-matched reads ( Figure 2). Almost 12% of the sRNAs are represented by the 24 nt class. Collectively, all mycelial libraries showed the majority of sRNA peaks at 24 followed by 26 nt, whereas in planta libraries showed the majority of sRNAs peaking at 26 followed by 24 nt.

Cluster analysis of mycelial and in planta libraries
We performed a cluster analysis on each M. oryzae supercontig, based upon the proximity-based algorithm of Lu et al. [30] for a comparative analysis of the mycelial and in planta small RNA libraries. This clustering algorithm allowed us to determine the abundance of sRNAs and their distribution across the M. oryzae genome by dividing the genome into a series of 500 base pair (bp) windows (clusters) and calculating the abundance of sRNAs for each cluster. The clusters are defined by fixed coordinates, which allows comparative analyses of multiple libraries and the identification of differentially regulated clusters (DRCs) in the genome. In our study, DRCs represent regions of the genome with a sRNA abundance that differs by 10-fold or more in pairwise comparisons between the libraries. All mycelial libraries were compared to CM, which resulted in a total of 477 DRCs. We examined the total number of sRNA clusters for each library and identified clusters that map to genic (which includes UTRs, introns and exons) and intergenic regions, various forms of repeat regions and transposable element regions. Out of the 477 total sRNA clusters, 251 were from genic regions, 213 were from intergenic regions and 20 were from repeat regions and transposable element regions ( Table 2). The CM vs. CS comparison revealed four DRCs, indicating carbon starvation had a minimal impact on the accumulation of sRNAs. NS vs. CM and PQ vs. CM had the greatest number of intergenic and genic-associated clusters. The genic and intergenic-associated clusters were 52 and 96 for NS vs. CM, and 179 and 69 for PQ   libraries; for the latter libraries, with the exception of two clusters mapping to retrotransposons, almost everything else associated with genic regions ( Table 2).
sRNA associate with different genomic regions, depending upon the fungal stress In total, taking data from all the libraries together, the majority of sRNAs were found to be associated with intergenic regions (54%), followed by repeat-containing regions (41%; Table 3). On the other hand, comparisons across libraries resulted in more DRCs associated with genes (Table 3). Although more total sRNAs mapped to intergenic regions compared to protein-coding regions, protein-coding regions contained a greater number of distinct sRNAs (Table 2). This difference may be explained by the presence of repeats which have not yet been annotated, and which were not identified by the RepeatMasker repeat identification program ( [31]; http://www.repeatmasker.org). Though approximately 50% of the genome codes for proteins and 7% is occupied by repeats [6], our sRNA data suggests only 4% of the total sRNAs sequenced matched to protein coding regions.
To determine whether our libraries identified repeatassociated siRNAs (known as rasiRNAs), we analyzed six different classes of repeats, and found that the majority of the sRNAs emanated from rRNA loci (89%) followed by retrotransposons (LTR/Gypsy and LINE/Tad1; 11%). Small RNAs that matched to transposons (DNA/TcMar-Pogo), simple repeats, and low complexity repeats occupied only a small fraction of the total sRNAs. Though rRNA-associated sRNAs are more abundant compared to retrotransposons-associated sRNAs, the distinct reads associated with the former were less (43%) compared to the latter (52%). The majority of sRNAs matching to protein-coding regions were in a sense orientation (77%; Table 3), and a similar phenomenon was observed for non-coding RNAs (tRNAs or rRNAs; 87%).
To further characterize the sRNAs associated with retrotransposons, representative loci from major repeat classes in M. oryzae such as LTR/Gypsy and LINE/Tad1 retrotransposons were analyzed for their size distribution, changes under different physiological conditions, strand specificity and 5' nt identity. Strand orientation of sRNAs was also included in our analysis because abundant antisense siRNAs are known to be associated with repeats [32] and are affected by the loss of RNAi machinery [5]. Also, studying strand specificity is important because sRNAs originating from overlapping gene pairs or natural cis-antisense gene pairs have been shown to play an important role in salt tolerance in Arabidopsis   [33] and gene regulation in mouse oocytes [34]. LTR/ Gypsy retrotransposons showed major peaks at 20-23 nt and abundant sRNAs in the antisense orientation (Additional file 2: Figure S2A and B). These sRNAs differentially accumulated under different physiological stress conditions; for example, one class of LTR/ GYPSY retrotransposon called GYMAG2 showed a spike in sRNA accumulation during nitrogen and oxidative stress conditions that was roughly 28% and 50% higher than any other condition (Additional file 2: Figure S2D). Analysis of the next most abundant major repeat class of retrotransposons, LINE/Tad1, revealed that its major peak for sRNAs based on abundance, was at 20 nt followed by 23/24 nt (Additional file 3: Figure S3A). Similar to the LTR/Gypsy retrotransposon class, sRNAs matching to LINE/Tad1 repeats also showed differential accumulation in various physiological stress conditions, but unlike the LTR class, LINE/Tad1 sRNAs did not show a spike under any one particular condition (Additional file 3: Figure S3C). Unlike LTR/Gypsy, sRNAs matching to LINE/Tad1 did not show abundant sRNAs in the antisense orientation (Additional file 3: Figure S3B). We observed that many of the annotated tRNAs and 5S rRNAs represent a significant portion of intergenic region -associated sRNAs. In order to study these genomic features, we selected random loci of each feature and analyzed their size distribution. The results showed a bimodal distribution for sRNAs originating from a majority of the rRNA loci, with size peaks at 20 and 23/ 24 nt (Additional file 4: Figure S4). Analysis of random tRNA loci revealed the presence of varied size distribution and increased sRNAs under the PQ condition (Additional file 4: Figure S4). We also observed strand bias for sRNAs originating from 5S rRNA and tRNA. All analyzed 5S rRNA and tRNA loci showed sRNAs on the sense strand (data not shown). sRNAs originating from mitochondrial sequences occupied only a small proportion, comprising 0.6% and 1.2% of total and distinct genome-matched reads, respectively (Table 3).
Stress conditions altered the proportion of individual classes of sRNAs associated with different genomic features (Additional file 5: Figure S5A). The PQ library revealed a notable reduction in both total and distinct abundance sRNAs matching to protein coding regions, compared to all the other libraries. Distinct abundance for the PQ library also showed an increase in sRNAs associated with intergenic regions, and with repeats. Another prominent difference was observed in the CS library; there was an increase in sRNAs associated with LTR Gypsy transposons, and a decrease in sRNAs associated with rRNA repeats (Additional file 5: Figure S5C). An increase in LTR retrotransposons was also observed during infection of rice plants with M. oryzae (Additional file 5: Figure S5D).

5' nucleotide preference
Nucleotide preference in the 5' position is the determining factor in plants for binding of sRNAs to specific AGO proteins [35]. While fungi do not appear to share this exact characteristic, preference for particular 5' nucleotides have been observed [12,21,36,37]. In M. oryzae, analysis of different libraries for their 5' nt preference revealed that sRNAs starting with U were abundant in both mycelial and in planta libraries, and there was a suppression for C in mycelial libraries and for G in in planta libraries (Figure 3). The 5' U preference was maintained in individual mycelial libraries -CM, MM, NS and PQ, but the CS library contained more sRNAs starting with G, closely followed by U. The suppression of C at the 5' position was also maintained in most of the mycelial libraries except for MM and PQ (data not shown). The 5' end nucleotide preference for the LTR class was U, and for LINE/Tad1 was both A and U (data not shown). The general pattern of abundant U and less G for sRNAs from in planta libraries was reflected in the LMg0 and LMg96, while LMg72 library had abundant A and less G (data not shown).

Correlation between transcriptional and posttranscriptional regulation under stress conditions
We hypothesized that transcriptional down-regulation of genes may correlate with a plethora of sRNAs accumulating either around that gene in the intergenic spaces, or within the gene itself. In order to gain a genome-wide perspective on how transcriptional and post-transcriptional gene regulation might be related, we compared our sRNA data to transcriptional genome data generated from M. oryzae under the same "stress" conditions as Mathioni et al. [29] in their microarray study. From this dataset, we identified 98 genes associated with clusters that showed differential accumulation of sRNAs (≥ 10-fold compared to CM control), 28 and 9 of which showed a negative correlation with microarray data under PQ and NS conditions, respectively (Additional file 6: Figure S6; Additional file 7: Table  S1; Additional file 8: Table S2; Figure 4). Further analysis of four of these genes showed that they have sRNAs originating from exons or introns in either the sense or antisense orientations (Table 4). Of these, only MGG_16711.7 (earlier annotated as MGG_05869.6) coding for a hypothetical protein containing an importin-beta N-terminal domain, had many sRNAs originating from its intron.
For the in planta libraries, we were most interested in sRNAs that were exclusively present at 72 hpi and/or 96 hpi (Additional file 9: Figure S7), indicating that these fungal sRNAs might be involved in invasive growth. In order to determine whether any of these sRNAs were negatively correlated with genes from the microarray study as described above, we first selected genes associated with sRNA clusters from the comparisons of LMg0 vs. LMg72, LMg0 vs. LMg96 and LMg72 vs. LMg96 that showed a ≥ 10-fold difference. Next, their corresponding loci from the CM vs. LMg72 data were compared to the microarray study [29]. This resulted in 12 genes that were both down-regulated in the microarray, associated with large quantities of sRNAs (Table 5). For each gene, abundant sRNAs present in LMg72 and/or LMg96 were analyzed for their position, and checked for their absence in non-inoculated rice sRNA libraries (http:// mpss.udel.edu/rice_sRNA), indicating that they were originating from the fungal, and not the rice, genome. The resultant sRNAs were largely found to be associated with exons, 5' UTR or 3' UTRs in a sense orientation (Table 5). Interestingly, one sRNA matching to MGG_08843 coding for the magnesium transporter, ALR2, was abundant at 72 hpi and present in an antisense orientation to the 5'UTR; Mathioni et al. [29] demonstrated down-regulation of this gene during carbon starvation, as well as in planta growth in rice and barley. Another gene includes the wellcharacterized transcription factor, ACE1 (MGG_04428.7); this gene is of interest, as it is known to be tightly regulated and involved in appressorial formation [38].
In order to confirm whether these genes might truly be under post-transcriptional regulation, we generated four sRNA biosynthetic mutants, and utilized real-time qRT-PCR to determine the expression of seven genes identified from the above analysis, namely MGG_0 8843.6, MGG_01596.6, MGG_04428.6, MGG_06609.6, MGG_05869.6, MGG_01439.6 and MGG_04470.6. Mutants tested included gene replacement deletions in the following loci: Dicer 1 and Dicer 2 (MGG_01541 and MGG_12357, previously generated and characterized by Kadotani et al. 2002, and hereafter referred to as MoDcl-1 and MoDcl-2), one RNA-dependent RNA Polymerase (RdRP MGG_13453; hereafter referred to as RdRP) and a double-knockout mutant in both Dicers (see Additional file 10: Figure S8 for deletion construct details, confirmation of deletions and Southern hybridization, revealing one disruption cassette per mutant line). We hypothesized that suppression of gene expression would be lifted if the genes for generating sRNAs in the fungus were no longer present, and sRNAs were no longer made. Six of the seven genes analyzed showed 1.5 to 1.9-fold increased expression over the wild type in the MoDcl1 mutant. The MGG_01596 gene, encoding a DNA damage response protein DUN1, increased in expression in the Δmodcl1, modcl2/modcl1 and in the Δmordrp-1 and was reduced in the Δmodcl2 mutant backgrounds. The MGG_04428 gene encoding the transcription factor ACE1 was increased~1.7-fold in the Δmodcl1 background ( Figure 5).

Discussion
Only a small fraction of RNA species (2-3%) are translated into proteins in higher organisms, while the remainder includes non-coding yet functional species such as small RNAs [39,40]. Small RNAs in eukaryotes have some well-established roles in modulating gene expression at the transcriptional and post-transcriptional levels, as well as in shaping the organization and modification of chromatin [41]. There is a robust body of knowledge on sRNAs in the model fungus Neurospora crassa, and in the budding yeast ( [5]; reviewed in [10,42]). Recently, several detailed studies on sRNAs have emerged in fungal pathogens, including M. oryzae amongst others [11,12,20,21,37]. The current study presented here, adds to this growing information by analyzing sRNAs under Our results demonstrate that the majority of the sRNAs originate from intergenic regions followed by repeats, protein coding regions and mitochondrial regions, when sRNAs from all libraries were taken together. Our results are supported by those of Nunes et al. [12], who reported abundant sRNAs matching to intergenic regions, transposable elements, rRNA, tRNA, and proteincoding genes, from mycelial and appressorial libraries.
Transposons account for about 7% of the M. oryzae genome (including 5.4% contribution by retrotransposon; [6]); our data suggests that only 4.5% of the total sRNA population matches to transposons (11% of the repeats population), which indicates RNAi may not act as the only mechanism to suppress transposons in M. oryzae. This observation is further supported by the hypothesis of Murata et al. [43], which supports the operation of a transposon-or repeat-silencing mechanism distinct from both transcriptional gene silencing and siRNA-mediated gene silencing.
After spore germination and appressorium formation, M. oryzae generates a penetration peg that enters into the plant's epidermal cell. Primary infection hyphae colonize this cell followed by formation of bulbous secondary hyphae and spreading to adjacent epidermal cells (the biotrophic phase; [44]). At about 72 hpi, secondary hyphae spread into mesophyll cells and at this time approximately 10% of total biomass is fungal [45], which could account for the reduced amount of total and distinct genome matched reads in in planta libraries, compared to the mycelial libraries. Under favorable disease conditions, visible lesions are formed by 96 hpi [45]. The drastic changes undergone by the pathogen during invasive, later stage growth are reflected in the increased number of genes showing differential abundance for sRNAs during 72 and 96 hpi, as determined by the cluster analysis (Table 2). During these later invasive stages, the fungus is believed to be coping with nutrient-limited environments [27,29]. In the in planta libraries, we observed 562 distinct reads associated with LMg72 compared to LMg96 after the removal of LMg0, which represent rice-only reads (Additional file 9: Figure  S7). Furthermore, out of the~8000-9000 sRNA reads found in both the LMg72 and LMg96 libraries, only 2132 overlapped, indicating that many unique small RNAs are generated during these invasive, later growth stages. We took advantage of previous transcriptional data generated from our lab in 2011 to determine   a Ratio of sRNA between rice library and CM library based on cluster data. b Expression of genes at 72hpi in rice compared to CM using microarray data [29]. c Abundance of signature.
whether clusters with sRNAs have a negative correlation with significantly differentially expressed genes in the microarray study. Ninety-eight such genes from the microarray study associated with sRNAs and 37.7% of these showed a negative correlation indicating that these particular genes may be regulated epigenetically. The remaining~60% of the microarray genes showed a positive correlation with the sRNA data. While these are relatively low numbers from a genome-wide perspective, the microarray data was generated from a very early version of the M. oryzae genome annotation, which has changed considerably in the last four years; furthermore, the 98 genes from the microarray study were based upon highly stringent criteria used both in this study, as well as the study from 2011, limiting the number of such genes we can identify.

Small RNAs differentially associate with retrotransposons, in a stress-dependent manner
Living organisms adapt to the environmental stress conditions by undergoing several phenotypic, physiologic and molecular changes, which requires precise regulation at the transcriptional, post-transcriptional and translational levels. These changes include but are not limited to activation and repression of several genes [46], signaling cascades [47], changes in molecular chaperones [48], activation of transposable elements [49,50] and chromatin changes and epigenetic modifications [51]. Several nutrient stress responsive sRNAs, have been identified in plants [52,53], bacteria [54] and Caenorhabditis elegans [55]. In pathogenic fungal species, efficient pathogen nutrition is essential for successful colonization and fungal fitness [56]. Transcriptome studies on M. oryzae and many other fungi have demonstrated that nutrient acquisition is a tightly regulated process, involving genes that function in amino acid metabolism and uptake, global nitrogen regulation and pathogenicity [27,[57][58][59]. In the present study, we found that repeatassociated sRNAs generally and LTR-retrotransposons specifically were induced by oxidative stress and nitrogen starvation, respectively (Additional file 5: Figure  S5). This observation may be similar to Arabidopsis roots undergoing phosphate starvation, which induces a high level of DICER-LIKE1-dependent sRNAs derived from the long terminal repeat of a retrotransposon and a group of 19 nt sRNAs that correspond to the 5' end of tRNA [60]. In rice, some repeat-derived sRNAs are down-regulated by abscisic acid and abiotic stresses accompanied by an increased accumulation of their precursors [61].
In M. oryzae, MAGGY retrotransposons are activated by heat shock and oxidative stress along with upregulation of MAGGY RNA [62] and an EST matching to a non-LTR-retrotransposon was identified in an M. oryzae-infected rice leaf library [63]. Similarly, our results revealed a spike in sRNAs associated with an LTR class during oxidative stress. Several ESTs matching to LTR-retrotransposons have also been identified in the in planta stage of the wheat fungal pathogen, Mycosphaerella graminicola, and are less abundant compared to their mycelial stage [64]. Similarly, our data revealed that distinct sRNAs reads matching to repeats in in planta growth were~10% higher than in the mycelial libraries (Additional file 5: Figure S5A and B). Reactivation of transposable elements during stress conditions have been found in many organisms from fission yeast to plants and animals (reviewed in [65]). Recently, a widely-studied tobacco transposon called TNT1 was found to be reactivated during infection of tobacco [66]. While we did see an increase in sRNAs associated with LTR/Gypsy retrotransposons during nitrogen starvation compared to the non-stressed condition, overall, there was not a striking difference in these types of associations in the case of the LINE/ Tad1 class of retrotransposon, where sRNA abundance was not markedly different in the stress versus non-stressed conditions. Together, our data and the wealth of information on transposon-associated small RNAs in other organisms suggest that retrotransposons might be suppressed via epigenetic mechanisms. Additional testing will be required to determine what underlies the increase in sRNAs associated with LTR elements during nitrogen starvation.
Interestingly, we noted an abundance of sRNAs on the antisense strands of LTR/Gypsy elements, while there was a roughly equal distribution of them on either strand of the LINE/Tad1 elements. This was different from the recent study by Nunes et al. [12] on sRNA species in M. oryzae, which showed a more equal distribution of siRNAs on both LTR-associated strands, and may be attributable to the differing conditions used in each study and differing sequencing depth. Association of small RNAs with one strand or the other has been wellstudied in Drosophila, for example, antisense siRNAs are abundant for repeat regions [32]. In general, transcription of transposons is activated by stress conditions; for example, carbon or nitrogen stresses activate the transcription of Drifter, a novel low copy transposon in Fusarium oxysporum. Some of the transposons are transcribed in both directions [67] and environmental stress-regulated antisense transcription of Ty1 LTRretrotransposons is observed in S. cerevisiae [68].
Derivation of the sRNAs in M. oryzae: antisense and 5' nucleotide position Around 25% of the tags that matched to genes are in antisense orientation (Table 3). Interestingly, one-fifth of the annotated genes of this fungus encode both sense and antisense transcripts [69]. The sources of the dsRNAs for Dicer processing include inverted repeats, bidirectional transcription and antisense transcripts from various loci [34]. There is a possibility that overlapping transcripts in an antisense orientation produced from a pair of neighboring genes on opposite DNA strands (convergent overlapping gene pairs or natural cis-antisense gene pairs) form dsRNAs, which may be further processed into sRNAs [33,34]. Recently, antisense transcripts were also found in other fungal species like the ectomycorrizal ascomycete Tuber melanosporum [70] and basidiomycete mushroom Schizophyllus commune [71]. It is also proven that siRNAs matching to mRNAs in S. castellii and in the RNAi reconstituted strain of S. cerevisiae were indeed produced from dsRNA formed by antisense and overlapping transcription [5,72]. Even in M. circinelloides, the majority of the exonic-siRNAS (ex-siRNAs) are antisense to mRNA [37]. Together, these results suggest that at least part of the siRNAs that matched to genes are produced from dsRNA formed by overlapping sense and antisense transcripts in M. oryzae. We observed that both 5S rRNA and tRNA loci were producing more of sense strand specific sRNAs whereas pseudo-tRNA loci have both sense and antisense sRNAs (data not shown). It has been shown that pseudo-tRNA of Bacillus cereus is detected to have a function outside protein synthesis [73]. Moreover, a subset of pseudogenes generates endo-siRNAs from double-stranded RNAs or directly in mouse oocytes [74]. The same scenario may exist in M. oryzae, because around 40% of pseudo-tRNAs are present either within a gene or opposite strand to it.
In the RNAi pathway, the Dicer-cleaved sRNAs are bound by Argonautes, and based on their 3' complementarity Argonautes exert homology dependent transcriptional cleavage or post-transcriptional silencing [75]. The binding of sRNAs to AGO is mainly dependent on structural features of precursor duplexes in animals [76,77] whereas in plants, it is determined by the 5' nt of the sRNA. In Arabidopsis AGO2 and AGO4 preferentially recruit sRNAs with 5' A, likewise AGO1 with U and AGO5 with C [35]. Abundant U at the 5' end observed in all the mycelial libraries in our study may be biologically significant (Figure 3), since this kind of preference was observed in other organisms including the single cell green alga Chlamydomonas [78]. Importantly, a preference for U (T) and suppression of C was observed in M. oryzae in a recent study, which our data supports [12]. In the fungal kingdom, 5' preference for U is observed in qde-2 (Argonaute) associated sRNAs and milRNAs, and disiRNAs of N. crassa [21]; ago-1 associated sRNAs of S. pombe [36]; sRNAs of S. castelleii and K. polysporus [5]; and Ex-siRNAs of Mucor circinelloides [37].
Though Dicer-like (DCL) proteins are evolutionarily conserved in eukaryotes, their numbers vary among organisms. In the multiple DCL protein-containing organisms such as Drosophila and A. thaliana, they function in distinct RNA silencing pathways with overlapping specificity [79]. Both M. oryzae and its close relative N. crassa encode two Dicers. However, MoDcl2 in M. oryzae functions in siRNA production under normal conditions, whereas N. crassa Dicers, orthologs of the M. oryzae Dicers, act redundantly. While it is likely that MoDcl1 does not have the ability to render or exert a silencing role under normal conditions either triggered by hairpin RNA-expressing transgenes or endogenous repetitive sequences, it is capable of producing siRNAs when overexpressed in an MoDcl2 deletion background. Moreover, its mRNA levels were up-regulated by 10-to 15-fold during the sexual stage. It has been proven that both MoDcl1 and MoDcl2 define functionally diversified RNAi pathways which may have arisen from transcriptional control and protein specialization [80].
In M. oryzae, deletion of MoDcl2 leads to increased abundance of MAGGY mRNA when it is introduced into the native genome [43]. With these reports in mind, it was expected that deletion of MoDcl2 in our experiments would increase the transcript abundance of genes which are known to possess more sRNAs in one or more libraries, in comparison with mycelia grown in complete media (control condition). Conversely, our results showed there is an overall reduction in transcript levels of these genes in the Δmodcl2 mutant and an increased level in the Δmodcl1 mutant. Interestingly, the known virulence gene ACE1 was among the few genes that fit this profile. ACE1 is tightly regulated at the transcriptional level; transcripts are only found in the penetration peg of the fungus, a specialized hypha that allows ingress through the surface of the leaf and into the plant cell [81]. We postulate that strong control of this gene might be governed at the posttranscriptional level. Further studies will reveal whether these and the other genes tested show the same profile when the biosynthetic mutants are subjected to the environmental stress conditions, or during pathogenicity. Several studies have shown that MoDcl2 is exclusively involved in siRNA production for transgenes and active transposons in M. oryzae, while less active transposons are not silenced by this same mechanism [43,80,82]. Moreover, MoDcl1 is also involved in silencing when it is overexpressed in a Δmodcl2 mutant background. Together, these findings along with results presented here suggest the hypothesis that MoDcl1, at least in M. oryzae, might define the existence of a MoDcl2-independent RNAi mechanism. The recent discovery of Argonaute associated Dicer independent siRNAs (disiRNAs) in N. crassa, the first discovered Dicer-independent siRNAs in the fungal kingdom, supports this hypothesis. Surprisingly, in M. circinelloides, Dicer2 plays a critical role in transgeneinduced silencing [25], but redundancy for the two Dicer genes was observed for production of endogenous sRNAs with DCL2 being the primary protein responsible for producing the majority of endogenous sRNAs [37]. Taken together, these studies suggest that M. oryzae presumably has similar redundant Dicer-based mechanisms for the production of sRNAs.

Conclusion
In conclusion, various stress-inducing and in planta conditions alter the sRNA profile of M. oryzae, with the main effect exerted on loci associated with protein coding regions. A subset of these genes showed an association with the transcriptome data, as well as a differential accumulation of their transcripts between the two Dicer-deficient mutant backgrounds. Stress conditions repressed the sRNAs associated with protein coding regions, and marginally increased the intergenic-and repeat-associated sRNAs. Considering our results and previous reports, questions of Dicer-independent mechanisms, and genetic redundancy of multiple copies of sRNA biosynthetic genes in the genome, warrant a detailed dissection of sRNA generating pathway(s) through the use of various sRNA biosynthetic mutants. Further studies on in planta interaction using sRNA biosynthetic mutants would provide better understanding of molecular mechanisms behind host-pathogen interaction which in turn useful for developing novel disease control strategies.

Methods
Fungal strain and growth conditions M. oryzae wild type strain 70-15 was used in this study because of the availability of genome sequence and transcriptome data. Filter paper stocks stored at −20°C were inoculated onto Oatmeal agar plates (OM: oatmeal 50 g/L and agar 15 g/L) for germination and to establish starter cultures. Fungus grown on OM plates were transferred to complete medium (CM: 10 g sucrose, 6 g yeast extract, 6 g casamino acid and 1 mL trace elements/L) plates and subsequently to complete liquid medium. Liquid grown fungus was harvested under sterile conditions and washed with sterile distilled water. The resulting fungal mat was carefully divided into equal pieces for further exposure to physiological conditions in liquid such as complete medium (CM), minimal medium (MM), carbon starved medium (CS: MM lacking carbon source), nitrogen starved medium (NS: MM lacking nitrogen source) and a reactive oxygen species-rich environment generated by amendment of Paraquat (PQ) to the CM, as described in detail by Mathioni et al. [29]. For the in planta libraries, twenty-day-old rice (Nipponbare) seedlings were mock inoculated with 0.2% gelatin and leaves were harvested immediately (LMg0). Conidia from M. oryzae 70-15 grown on OM plates for 10 days (d) were harvested using 0.2% gelatin. Concentration of the conidial suspension was adjusted to 1 × 10 5 spores/mL and 3 mL of the suspension was spray inoculated onto rice plants; leaves were harvested after 72 hours post-inoculation (hpi) and 96 hpi respectively for LMg72 and LMg96 libraries.

Small RNA library construction and data analysis
Total RNA was isolated using Trizol reagent (Invitrogen, New York, USA) according to manufacturer's instructions. The quality of the RNA was checked by agarose gel electrophoresis and NanoDrop ND-1000. High molecular weight RNA was removed by PEG precipitation, and low molecular weight RNA was size-fractionated on a polyacrylamide gel. Gel portions corresponding to 20-30 nt were recovered and small RNA libraries were constructed by following the protocol by Lu et al. [83] with modifications using standard Illumina sRNA adapters. The quality of the resulting libraries was checked by subjecting an aliquot of each library to TOPO® cloning (Invitrogen, USA), and approximately 96 clones from each library were sequenced by the Sanger method. The resulting sequences were checked for the presence of adapters, size of cDNAs, and their genome match. The libraries were then sequenced using Illumina's Genome Analyzer II (GA II ) at the Delaware Biotechnology Institute. Adapters were removed from the resulting sequences and matched to the M. oryzae genome (version 6 assembly) from the Broad Institute, and reads with sequencing errors were filtered out using PERL scripts. tRNA, rRNA, snRNA or snoRNAsmatched reads were removed from the total reads, and genome-matched reads were obtained using the program Bowtie, accepting only perfect matches. From the total genome matched reads, distinct genome matched reads were calculated. The abundance of each read was normalized to 5 million units of transcripts per 5 million (TP5M) to account for the difference in the sequencing depth and to facilitate accurate comparisons across libraries. To compare the different libraries, the entire genome was divided arbitrarily into 500 bp windows/clusters based on the proximity-based algorithm ( [30]; described in the Results section). Clusters showing ≥ 10-fold difference between 2 or more libraries were selected for further cluster analysis. Repeat analysis was performed using programs namely, RepeatMasker [31], etandem and einverted [84].

Targeted disruption of genes involved in sRNA biogenesis and mutant screening
Gene replacement cassettes were constructed using adaptamer-mediated PCR [85]. Typically, 1.3 kb of upstream and downstream sequence of each target gene was amplified with primers that contained adaptamer sequences. A 1.4 kb fragment containing the hygromycin B phosphotransferese gene with the trpC promoter from A. nidulans was amplified from plasmid pCB1003 using the adaptamer sequence attached to the forward HPHF and reverse HPHR primer set (Additional file 11: Table  S3). Using nested primers, three individually amplified fragments of upstream and downstream sequence and hygromycin resistance gene were combined and amplified together to construct a hygromycin cassette for gene replacement of approximately 3.3 kb in length. The hygromycin cassette was transformed into 70-15 protoplasts as previously described [86]. Gene disruption mutants were identified by PCR screening using primers outside the flanking regions and gene specific primers (Additional file 11: Table S3) and further confirmed by Southern blot analysis (Additional file 10: Figure S8). For deletion of MoDcl1 and MoDcl2, gene replacement cassettes were constructed as described previously with the modification of BAR gene, amplified from pCB1530, being used as the selectable marker [87] in place of HYG. Gene replacement cassettes for MoDcl1 and MoDcl2 were transformed into protoplasts of MoDcl2 and MoDcl1 knockouts, respectively.
Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) Genomic DNA was removed from total RNA by TURBO® DNase treatment (Ambion, New York, USA). Reverse transcription reactions were performed with 500 ng of RNA in 20 μl reaction using High Capacity RNA-to-cDNA Kit (Applied Biosystems, California, USA) by following the manufacturer's instructions. Real time PCR was done using 5 PRIME SYBR mastermix for Real-time PCR (5 PRIME, Maryland, USA) in an Eppendorf Mastercycler using 2 μl of diluted reverse transcribed templates. All samples were repeated three times. Relative expression values were calculated using the 2 -ΔΔCT method with GADPH (MGG_01084.6) as a housekeeping control. Primers used are listed in Additional file 12: Table S4.

GEO accession numbers
The data obtained from small RNA deep sequencing studies were deposited in the Gene Expression Omnibus database at NCBI (http://www.ncbi.nlm.nih.gov/geo/). The accession number for GEO is GSE43277). The data can also be found at http://mpss.udel.edu/mg_sRNA.