Transcriptional and Small RNA Responses of the White Mold Fungus Sclerotinia sclerotiorum to Infection by a Virulence-Attenuating Hypovirus.

Mycoviruses belonging to the family Hypoviridae cause persistent infection of many different host fungi. We previously determined that the white mold fungus, Sclerotinia sclerotiorum, infected with Sclerotinia sclerotiorum hypovirus 2-L (SsHV2-L) exhibits reduced virulence, delayed/reduced sclerotial formation, and enhanced production of aerial mycelia. To gain better insight into the cellular basis for these changes, we characterized changes in mRNA and small RNA (sRNA) accumulation in S. sclerotiorum to infection by SsHV2-L. A total of 958 mRNAs and 835 sRNA-producing loci were altered after infection by SsHV2-L, among which >100 mRNAs were predicted to encode proteins involved in the metabolism and trafficking of carbohydrates and lipids. Both S. sclerotiorum endogenous and virus-derived sRNAs were predominantly 22 nt in length suggesting one dicer-like enzyme cleaves both. Novel classes of endogenous small RNAs were predicted, including phasiRNAs and tRNA-derived small RNAs. Moreover, S. sclerotiorum phasiRNAs, which were derived from noncoding RNAs and have the potential to regulate mRNA abundance in trans, showed differential accumulation due to virus infection. tRNA fragments did not accumulate differentially after hypovirus infection. Hence, in-depth analysis showed that infection of S. sclerotiorum by a hypovirulence-inducing hypovirus produced selective, large-scale reprogramming of mRNA and sRNA production.

1. Introduction

Fungal viruses (mycoviruses) are highly diverse, infect pathogenic and nonpathogenic fungi, and can significantly reduce the virulence of pathogenic fungi [1]. Fungal–virus interactions involve the interplay between gene expression networks, some of which are influenced by small noncoding RNAs [2]. Central to those interactions are RNA-mediated antiviral defenses that are activated by double-stranded RNA (dsRNA) in a process termed RNA silencing [3,4,5]. Specificity of the defense is imparted by short (21–24 nt) interfering RNAs (siRNAs) produced from viral dsRNA by RNase III-type enzymes called dicer-like (DCL) proteins [6]. One strand of each siRNA is combined with Argonaute (AGO) proteins into RNA-induced silencing complexes and directs degradation of complementary viral RNA sequences [3].

While RNA silencing genes are highly conserved in animals and plants, they are much less conserved in fungi and oomycetes [7]. For example, the genomes of Saccharomyces cerevisiae and Ustilago maydis lack homologs of canonical AGO and DCL genes, while they are present in the genomes of Saccharomyces castellii and Ustilago hordei [8]. In addition, some fungi, e.g., Candida albicans, have evolved novel DCL proteins that differ significantly from those of higher eukaryotes [9], making it challenging to ascertain the importance of fungal endogenous small RNA processing in response to virus infection.

Endogenous small RNAs (sRNAs) play important roles in gene regulation in development and responses to biotic and abiotic stresses in animals and plants [10], but their roles are not well defined in phytopathogenic fungi. It remains elusive whether filamentous fungi have functional microRNAs (miRNAs), typically 21–24 nt in length, generated by DCL processing of short hairpin structures and regulate gene expression directly through interactions with AGO and indirectly through the production of secondary phased siRNAs (phasiRNAs) from long noncoding RNAs [11,12,13]. PhasiRNAs are a class of secondary sRNAs triggered by miRNAs; the substrate of phasiRNAs is produced by host RNA-dependent RNA polymerase (RDR). The presence of phasiRNAs is an indication of functional miRNAs.

While most plants share a common set of miRNAs [10], fungal species analyzed so far however express diverse suites of miRNAs [14,15]. A few studies have attempted to identify the roles fungal miRNAs play in regulating gene expression. Small RNA and high-throughput rapid amplification of cDNA ends (HT-RACE) to identify cleavage targets for Fusarium oxysporum microRNA-like RNAs (milRNAs) did not identify any milRNAs that were present in the current miRNA database, and none of the milRNAs were predicted to trigger cleavage of F. oxysporum mRNAs [16]. Analysis of differentially expressed sRNA loci and miRNA accumulation in Aspergillus flavus under different growth conditions suggested that miRNAs play important roles in cellular functions including mycotoxin biosynthesis and mycelial growth [17]. The predicted targets for sRNAs from the wheat stripe rust fungus (Puccinia striiformis) were enriched for kinases and small secreted proteins suggesting that development-related signaling pathways are regulated by sRNAs in P. striiformis [14].

Another source of abundant sRNAs is derived from tRNAs. In plants and animals, mature tRNAs have been confirmed to be sources of functional sRNAs produced in a DCL-independent manner that have been implicated in post-transcriptional and epigenetic regulation of gene expression and repression of retrotransposons [18,19,20]. In the plant-pathogenic fungus Magnaporthe oryzae, tRNA-derived RNA fragments (tRFs) were more abundant in appressoria—specialized-infection tissues—than in vegetative mycelia [21], suggesting that small RNAs in fungi also play active roles in the regulation of growth and development as in higher eukaryotes.

Virus-derived siRNAs (vsiRNAs) have been shown to direct the cleavage of host RNAs in plants [22,23,24,25]. Also, another new class of endogenous small RNAs were shown to be activated in Arabidopsis plants in response to infection by cucumber mosaic virus that were termed virus-activated siRNAs (vasiRNAs) that lead to the silencing of a broad set of host genes and establishment of an antiviral state [26]. Clearly, small RNA are not just products of RNA silencing but instead exhibit biological functions to down regulate gene expression.

Changes in fungal gene expression in response to virus infection vary greatly depending on the nature of the host–virus interaction that hinges on the activity of antiviral RNA silencing pathway. Saccharomyces cerevisiae LA virus, a member of the Totiviridae, and its M1 satellite impart a beneficial killer phenotype to S. cerevisiae, which lacks RNA silencing genes, and produces relatively small changes in gene expression [27]. In contrast, infection of plant pathogenic fungi with robust RNA silencing by mycoviruses that alter fungal virulence have much more pronounced effects on fungal gene expression. For example, infection of Cryphonectria parasitica with Cryphonectria hypovirus 1 (CHV1) significantly altered the expression of more than 13% of the analyzed transcripts [28]. Proteomic analysis of proteins secreted by CHV1-infected C. parasitica identified 99 proteins with differential accumulation relative to virus-free cultures [29]. Similar proteomic analysis of infection of Fusarium graminearum by Fusarium graminearum virus DK21, which perturbs development and attenuates the virulence of the fungus, altered the accumulation of nearly 150 proteins [30].

High-throughput (HT) sequencing of different RNA species provides an approach to functionally characterize how small RNAs interact with potential targets to determine whether milRNAs and siRNAs repress gene expression by mRNA cleavage in fungi. The combination of three types of sequencing data—small RNA-seq, HT-5′-RACE (degradome), and mRNA-seq—allows a global analysis of small RNA function. By analyzing predicted milRNA sequences together with mapped degradome reads, high confidence cleavage events can be identified. Predicted cleavage events can be supported by observations of downregulation of predicted mRNA targets in RNA-seq analysis. Confirmation of predictions is especially important in fungi because, unlike plant or animal miRNAs, fungal pre-milRNAs are do not have well-defined secondary structures and the thermodynamics of the interactions between miRNAs and targets is not well defined, which can lead to false positive cleavage site predictions that cannot be confirmed empirically.

Because mycoviruses can reduce fungal virulence to animal and plant hosts, we aimed to examine the effect of mycovirus infection on gene expression and different classes of sRNA accumulation and gather functional evidence for the existence of miRNAs in the white mold fungus Sclerotinia sclerotiorum, a devastating plant pathogen. Previously, we determined that Sclerotinia sclerotiorum hypovirus 2-L (SsHV2-L) reduced/delayed development in sclerotia, enhanced production of aerial mycelia, and induced hypovirulence. Now, we report the use of high throughput sequencing to (1) show that infection by SsHV2-L changes in mRNA and sRNA accumulation in S. sclerotiorum; (2) identify new classes of fungal sRNAs; and (3) examine whether predicted sRNA cleavage events correspond to changes in gene expression. Supported by HT-RACE analysis, we show that endogenous S. sclerotiorum sRNAs were capable of directing the cleavage of coding and noncoding RNAs, the latter leading to the production of phasiRNAs, demonstrating the evidence of miRNA activity.

2. Materials and Methods

2.1. Preparation of Sclerotinia sclerotiorum Cultures and RNA Extraction

Virus-free (VF) and virus-transfected (VT) cultures of S. sclerotiorum strain SsDK3 were produced as described previously [31]. Two sequencing runs were performed: Trial #1 and Trial #2. For the sequencing run of 4 libraries (Trial #1), fungal cultures were grown in potato dextrose broth (PDB) at 25 °C for 10 d. For the separate sequencing run of 10 libraries (Trial #2), fungal cultures were grown on potato dextrose agar (PDA) at 21 °C for 4 d. Total RNAs were extracted from VF and VT cultures using RNeasy Mini (Qiagen, Valencia, CA, USA). Small RNAs were extracted from two biological replications each for VF and VT cultures using mirVana miRNA isolation kits (ThermoFisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions.

2.2. Analysis of S. sclerotiorum Transcriptome

Trial #1 of four RNA-Seq libraries (two VF and two VT) was prepared using the TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA). The RNA-Seq libraries were barcoded and sequenced as single-end 100-nt reads in one lane on an Illumina HiSeq2500. Trial #2 of ten RNA-Seq libraries (five VF and five VT) were prepared using the same kit and barcoded/sequenced as single-end 100-nt reads on an Illumina HiSeq4000. All RNA-Seq reads were assembled de novo using Trinity [32] and used to identify novel coding regions in the S. sclerotiorum genome using BRAKER1 [33]. The abundance of reads aligning to predicted coding regions was estimated using RSEM [34]. Differentially expressed coding regions were identified using the DESeq2 Bioconductor package [35]. Differentially expressed coding regions were classified using PANTHER [36]. The relative abundance of transcripts from S. sclerotiorum coding regions was determined using reads per kilobase of coding sequence per million mapped reads (RPKM) [37].

2.3. Analysis of S. sclerotiorum Small RNA Populations

Four sRNA libraries (two VF and two VT) were prepared using the TruSeq Small RNA Library Prep Kit (Illumina). The four sRNA libraries were barcoded and sequenced in one lane as 50-nt single-end reads on an Illumina HiSeq2500. Loci with significantly different accumulations of sRNAs between VF and VT samples were identified using DESeq2. Candidate miRNAs were predicted with six different programs, miRDeep2 [38], MiRDeep* [39], miRDeep-P [40], miReap (http://sourceforge.net/projects/mireap/), MiRPlant [41], and ShortStack [42], using the default settings for each program. Small RNA sequences with fewer than 20 reads, and sequences that aligned to coding regions or ribosomal RNA or more than 10 locations in the S. sclerotiorum genome were excluded. Loci producing phased sRNAs were identified with Shortstack [43]. Plots of phased siRNAs were prepared using the equation described by Howell et al. [44]. Candidate tRNA genes in the S. sclerotiorum genome were predicted using the tRNAscan-SE server [45].

2.4. High-Throughput RNA Ligase-Mediated Rapid Amplification of cDNA Ends (HT-RACE) Analysis

Four HT-RACE (degradome) libraries (two VF and two VT) were constructed as described by Li et al. [46]. Briefly, the Illumina sRNA-seq 5′ adapter was ligated to the purified polyadenylated RNAs from two VF and two VT cultures using T4 RNA ligase. The 5′ termini of most intact fungal mRNAs are blocked by an m7GTP cap structure, while cleaved RNAs have a 5′ terminal phosphate that is available for ligation with the adapter. The 5′-adapter-ligated RNAs were purified, fragmented, treated with phosphatase, purified, and ligated to the Illumina sRNA-seq 3′ adapter. The dual-adapter ligated RNAs were reverse transcribed using the Illumina sRNA-seq RT primer and amplified with the Illumina Gx1 and GX2 sRNA-seq PCR primers. The four barcoded libraries were pooled and sequenced in one lane for 50-nt single-end reads. Potential targets for candidate miRNAs and cluster sequences were identified using the HT-RACE sequence data and the predicted S. sclerotiorum coding sequences with the CleaveLand4 pipeline [47]. Predicted miRNA targets with P-values of less than 0.05 and degradome categories 0 and 1 hits were selected for downstream cluster analysis.

3. Results

3.1. Changes in mRNA Accumulation Associated with Infection of S. sclerotiorum by SsHV2-L

Two sequencing runs of the transcriptome were performed. Trial #1 had two biological replications each, and Trial #2 had five biological replications each. The accession numbers for the infected and virus free treatments from the Trial #1 are SRR8306347 and SRR8306348, respectively, while SRR8305679 and SRR8305680 are the accession numbers for the infected and virus free treatments from the Trial #2, respectively.

RNA-Seq analysis of Trial #1 produced a total of 2.11 × 108 100-nt reads from the four libraries. On average, 89.9% and 87.1% of the reads from the virus-free and SsHV2-L-infected libraries aligned to the S. sclerotiorum genome sequence [48], respectively (Table S1a), and assembled into 23,968 contigs with an N50 of 2834 nt. RNA-Seq analysis of Trial #2 produced a total of 3.99 × 108 100-nt reads from the ten libraries with 89.0% and 83.7% of the reads from the virus free and SsHV2-L infected libraries aligned to the S. sclerotiorum genome sequence, respectively (Table S1b). The S. sclerotiorum genome is predicted to contain 14,053 coding regions [48], but re-annotation of the S. sclerotiorum genome using the RNA-Seq data identified an additional 174 putative coding regions with an average size of 680 nt in this study. Overall, in the RNAs from the infected samples, an average of 1.45% and 5.7% of the sequence reads were derived from the SsHV2-L genome from Trial #1 and Trial #2, respectively.

The RNA-Seq data analyzed from Trial #1 when the cultures are 10-d old in PDB at 25 °C, among the proteins differentially expressed primarily involved in carbohydrate metabolism, Sclerotinia sclerotiorum hexose transporter 1 (Sshxt1; SS1G_04841), and a predicted invertase (SS1G_07184; beta-fructofuranosidase) were upregulated 37.8-fold and 18.8-fold, respectively, in SsHV2-L-infected cultures. The expression of S. sclerotiorum hexose transporter 2 (SS1G_13734) was upregulated 89.3-fold by SsHV2-L infection. In Trial #1, 958 coding regions differed significantly; 471 were upregulated and 487 were downregulated (Table S2a). The same trial also found that SsHV2-L infection significantly enhanced the accumulation of mRNAs predicted to encode a DEAD-box ATP-dependent RNA helicase (SS1G_04322) and a mitochondria-localized pentatricopeptide repeat-containing protein (SS1G_08638) (Table S2a). Similar to plant gene expression changes induced by plant virus infection [49], mycovirus SsHV2-L infection also increased transcription levels of retrotransposon-related coding sequences in the fungus. SsHV2-L enhanced expression of coding regions for ubiquitin-related proteins (SS1G_02395 and SS1G_00267), aldo/keto reductase (SS1G_00727), and a heat shock protein (SS1G_12897). SsHV2-L infection altered the accumulation of mRNAs encoding 10 cytochrome P450s and 13 methyltransferases. SsHV2-L downregulated SS1G_01530 which is a putative SPX-domain protein in a manner similar to the S. cerevisiae protein SPX-domain protein Syg1 that interacts with the β subunit of G-proteins to suppress the lethality of α subunit deficiency [50].

The two trials of transcriptome sequencing runs from S. sclerotiorum grown on different media, for different times and at different temperatures generated very different lists of differentially expressed genes in response to SsHV2-L infection. Both sequencing runs showed that infection by SsHV2-L did not significantly alter the expression of S. sclerotiorum homologs of AGO (SS1G_00334 and SS1G_11723), DCL (SS1G_13747 and SS1G_10369), or RDR (SS1G_03377, SS1G_13161 and SS1G_09915). With the ten libraries of Trial #2 RNA-Seq analysis from 4-d old cultures grown on PDA at 21 °C, infection of S. sclerotiorum with SsHV2-L significantly (p < 0.05) altered the abundance of transcripts from 1319 coding regions; 779 were upregulated and 540 were downregulated. Panther analysis of the functional classifications of the annotated S. sclerotiorum coding regions found that sequences predicted to be involved in rRNA metabolic processes (P = 6.4 × 10−20) and cellular component biogenesis (P = 4.6 × 10−13) were significantly over represented in the SsHV2-L-infected cultures. rRNA-processing protein 7 homolog A-related protein (SS1G_00995) was downregulated by 2.8-fold. Ribosomal U3 small nucleolar ribonucleoprotein proteins IMP3—SS1G_05011 and SS1G_05012—were both downregulated ~2.5-fold. SS1G_03709 and SS1G_00849, which are predicted to be virulence effectors when infecting soybean [51], were downregulated by 5.6-fold upon SsHV2-L infection. Several methyltransferases were downregulated, including SS1G_12790 by 4.9-fold and SS1G_11246 and SS1G_01682 both by 4-fold (Table S2b). Inactive methyltransferase may lead to terminal modification of small RNAs that have been found to trigger phasiRNA production in plants [52].

3.2. Small RNA Accumulation in Healthy and SsHV2-L-infected S. sclerotiorum Cultures

Sequence analysis of the four sRNA libraries produced a total of 9.04 × 107 reads of 17 to 36 nt, of which, 63.5% aligned to the S. sclerotiorum genome sequence (Table S3). The accession numbers are SRR8306349 and SRR8306350 for the virus-infected and virus-free treatments, respectively. On average, 16.5% of the sRNA reads aligned to the 192 S. sclerotiorum tRNAs predicted by tRNAscan-SE [45], 12.9% to ribosomal RNAs, 5.9% to retrotransposons, and similar to Fusarium oxysporum and N. crassa, 12.9% of the S. sclerotiorum sRNA reads aligned to the S. sclerotiorum mitochondrial genome [16,53,54]. In the RNAs from the two infected samples 9.0% of the sequence reads were derived from the SsHV2-L genome on average.

In all four samples, 22-nt sRNAs were the most abundant with a strong preference for uracil at the 5′ position (Figure 1B). Among the reads that aligned to unique positions in the S. sclerotiorum genome, there was a second peak at 27-nt (Figure 1A). Most of the 27-nt sRNAs aligned to S. sclerotiorum ribosomal RNA genes (Figure 1D). Few of the 27-nt sRNAs aligned to coding regions or sequences with homology to retrotransposons (Figure 1C). The density with which reads aligned to the genome sequence was highly variable. Twenty of the 25 most abundant nonribosomal RNA and noncoding sequences were of 27- to 32-nt in length. Most aligned to three to six loci in the S. sclerotiorum genome and many were conserved in the genome of B. cinerea. For example, the most abundant sequence (5′-UCCGAAUUAGUGUAGGGGUUAACAUAACUC-3′) with over 4.5 × 106 reads (5.0% of total reads) aligned to one locus each on S. sclerotiorum chromosomes 4 and 7, and two loci on chromosome 16 that corresponded to the locations of predicted glutamine tRNAs (Figure 1E). The sequence, designated as tRF5-Glu(GAA), also aligned to predicted glutamine tRNAs on B. cinerea chromosomes 2, 8, and 14 and two loci on chromosome 16 (data not shown). Hence the highly abundant sRNAs were likely derived from mature tRNAs and probably represented tRNA halves that are produced by endonucleolytic cleavage of mature tRNAs, sometimes in response to biotic and abiotic stresses [18]. However, the abundance of sRNAs derived from tRNAs was on average 1.4-fold lower in infected-SsHV2-L samples than virus-free controls.

Figure 1

Size distributions of small RNA sequences that aligned to the S. sclerotiorum genome. (A) Size distribution of small RNA libraries from combined virus-free and hypovirus-infected S. sclerotiorum cultures. White and black columns represent unique and total reads of the sRNAs, respectively. (B) Frequency of 5′ terminal nucleotides from pooled small RNA samples. Size distribution of small RNA reads aligning to (C) coding regions, intergenic regions, retrotransposon sequences, and (D) ribosomal RNA and tRNA sequences. (E) Mature tRNA structures predicted by tRNAscan-SE with sequences of the two most abundant small RNA sequences that resembled stress-induced tRNA halves. I: 4.5 × 106 reads (5.0% of total reads); tRNA Glu-derived tRF5-Glu(GAA) on S. sclerotiorum chromosomes 4, 7, and 16 (two copies); and B. cinerea chromosomes 2, 8, 14, and 16 (two copies). II: 1.6 × 106 reads (1.9%); tRNA Asp on S. sclerotiorum chromosomes 1, 5, 11, 12, and 14; and B. cinerea chromosomes 9, 10, and 13 (two copies).

Similar to sRNA reads derived from the S. sclerotiorum genome, the predominant size of sRNAs derived from the SsHV2-L genome were 22 nt in length (Figure 2A) and nearly 90% of the aligned reads contained a 5′-terminal U residue (Figure 2B). As has been reported for other virus infections, the sRNA reads were derived nonuniformly from both RNA strands (Figure 2C).

Figure 2

Size distribution of small RNA sequences that aligned to the Sclerotinia sclerotiorum hypovirus 2 L (SsHV2L) genome sequence. (A) Size distribution and (B) frequency of 5′ terminal nucleotides of small RNAs that aligned to the SsHV2L genome. (C) Distribution of small RNA reads that aligned to the SsHV2L genome. Bars above zero indicate alignment to the positive strand, and bars below zero indicate alignment to the negative strand.

3.3. Identification of S. sclerotiorum Loci Producing microRNA-Like RNAs

Cluster analysis in ShortStack identified 8617 intergenic loci that produced at least 100 sRNA reads. Nearly 10% (835) of the intergenic loci producing sRNAs responded significantly to SsHV2-L infection. Because of the variability in RNA silencing-related genes in fungi, candidate miRNAs were predicted from the sRNA data using six different programs. The analysis predicted a total of 459 candidate miRNAs in noncoding and nonribosomal regions of the S. sclerotiorum genome with at least 20 reads in the combined data set. Notably, the numbers of candidate miRNAs predicted differed for each of the ten programs used because of the evolution position of fungi between plants and animals. For example, miRDeep2 [38], MiRDeep* [39], miRDeep-P [40], miReap (http://sourceforge.net/projects/mireap/), MiRPlant [41], and ShortStack [42], predicted 46, 173, 175, 27, 139, and 0 candidate miRNA, respectively (Table S4). Just 18% of the candidate miRNAs were predicted by more than one of the programs, among which the predicted structures of candidate miRNAs 0159, 0714, and 1287 contained asymmetrical bulges (Figure 3A).

Figure 3

Examples of sRNA producing loci in the S. sclerotiorum genome (A) capable of folding into structures similar to pre-microRNA and (B) conserved in the genomes of other members of the family Sclerotiniaceae. Arrows indicate the positions of mature microRNA-like sequences. Connected shaded boxes indicate regions of conserved base pairing in predicted stem-and-loop structures.

Because animals, plants, and fungi produce lineage-specific and species-specific miRNAs [10,55,56], the conservation of the candidate miRNA loci of S. sclerotiorum was analyzed in the family Sclerotiniaceae. Among the candidate miRNA loci, just 18.6% were conserved in at least one other member of the family Sclerotiniaceae. For example, candidate miRNAs 0386 and 0437 were conserved in the genomes of Botrytis cinerea and Sclerotinia homoeocarpa (Figure 3B).

3.4. Potential Targets for S. sclerotiorum sRNAs

Because in silico prediction of miRNA targets alone can produce a large number of false positives in fungi [57], HT-RACE or degradome was used to identify mRNA with uncapped 5′ termini that could represent cleavage events. The abundance of HT-RACE reads were mapped on S. sclerotiorum coding sequences relative to the positions of peak cluster sequences and candidate miRNAs sequences with CleaveLand [47] to detect possible sRNA-mediated slicing sites. CleaveLand separated potential cleavage sites into five categories depending on the number of HT-RACE reads that aligned relative to the sRNA. In category 0 cleavage sites, the sRNA aligned over the most abundantly detected 5′ terminus in the coding sequence. Eighteen candidate miRNAs and 22 sequences from cluster analysis had degradome P values of less than 0.05 and category 0 or category 1 slicing sites (Table S5). As has been reported in other systems [58], one of the milRNAs, which was also identified by the cluster analysis, was predicted to direct cleavage of transcripts of SS1G_00334, which is predicted to encode an Argonaute 2 homolog. Other predicted targets included an ammonium transporter (SS1G_04502), an extracellular membrane protein (SS1G_12056), a pachytene checkpoint protein (SS1G_12812) and a retrovirus-related polymerase (SS1G_12103) (Table S5).

3.5. Detection of Phased siRNAs from S. sclerotiorum Noncoding RNAs

Most (78.8%) of the S. sclerotiorum sRNAs were derived from intergenic regions, and de novo assembly of the S. sclerotiorum transcriptome data identified 8761 noncoding transcripts. Analysis of the distribution of phasing of sRNA reads in the S. sclerotiorum genome identified 164 loci producing 21-nt or 22-nt phasiRNAs with phasing scores greater than 25 (Figure 4A). Thirty-six loci showed significantly different (p < 0.05) accumulations of sRNAs between mock-inoculated and SsHV2-L-infected samples (Figure 4A). Seventeen showed reduced accumulation and 19 showed enhanced accumulation of sRNA sequences. Plots of phasing scores of sliding 10-cycle windows showed characteristic repeating patterns for phased loci (Figure 4B). The detection of large numbers of loci producing phasiRNAs suggests that S. sclerotiorum produces sRNAs that function similar to miRNAs in higher eukaryotes.

Figure 4

Genome-wide identification of loci producing phased small RNAs in S. sclerotiorum. (A) Distribution of loci producing phased small RNAs (red horizontal bars) on S. sclerotiorum chromosomes. Loci marked with asterisks showed significantly different accumulation of small RNAs between mock-inoculated and SsHV2L-infected samples. Asterisks on the left of the bars indicate reduced small RNA accumulation; bars to the right indicate increased small RNA accumulation. (B) Small RNA abundances and phasing score distributions across two loci producing phased small RNAs in the S. sclerotiorum genome.

4. Discussion

In this study, we analyzed changes in mRNA and sRNA accumulation in S. sclerotiorum in response to persistent hypovirus infection and found that infection by SsHV2-L altered the accumulation of nearly 10% of coding mRNAs and sRNAs. The ability to experimentally inoculate S. sclerotiorum with SsHV2-L using in vitro transcripts allowed us to fine map virus-induced changes in mRNA and sRNA accumulation during hypovirus infection. Related studies examined changes in gene expression induced by hypovirus infection in C. parasitica and F. graminearum. In C. parasitica, CHV1 infection altered the expression of mRNAs predicted to encode proteins involved in carbon metabolism, stress responses, and regulation of transcription [28], which was supported by metabolomics analyses [59]. Even though FgHV1 does not induce hypovirulence, FgHV1 infection of F. graminearum also significantly altered the accumulation of mRNAs predicted to encode proteins involved in sugar and carbohydrate metabolism [60]. Metabolism of mannitol and other sugars are important in compatible fungal–host interactions [61]. Hence, persistent SsHV2-L infections may be metabolically similar to latent viral infections where carbon is directed away from the TCA cycle for glycolysis to reduce apoptosis [62].

Some of the transcripts differentially expressed in this study provided insights on evolutionarily conserved basal defense systems. Infection of S. sclerotiorum by SsHV2-L induced expression of mRNAs predicted to encode stress-responsive proteins including a putative mitochondria-localized pentatricopeptide repeat-containing protein and a DEAD-box RNA helicase. In Arabidopsis, the stress-induced mitochondria-localized pentatricopeptide repeat protein, PGN, is involved in defense against necrotrophic fungi and tolerance to abiotic stresses [63]. In addition to their deeply conserved roles in stress responses [64], host DEAD-box RNA helicases enhance and sometimes are required for the replication of animal and plant RNA viruses [65,66]. Hence, it is possible that S. sclerotiorum DEAD-box RNA helicases could enhance SsHV2-L accumulation.

Infection of S. sclerotiorum by SsHV2-L did not alter the accumulation of mRNAs encoding RNA silencing-related proteins. In contrast, infection of C. parasitica by CHV1 upregulated DCL2 mRNA accumulation up to 20-fold [67]. One of the four C. parasitica AGO genes, AGO2, which is required for antiviral defense, was also upregulated by CHV1 infection [68,69]. Similarly, one DCL, one AGO, and two RDR genes were upregulated in Rosellinia necatrix by infection with Rosellinia necatrix mycoreovirus 3 or Rosellinia necatrix megabirnavirus 1 [70]. Similar to infection of S. sclerotiorum by SsHV2-L, infection of R. necatrix by Rosellinia necatrix partitivirus 1, Rosellinia necatrix quadrivirus 1, and Rosellinia necatrix victorivirus 1 [70] did not significantly alter the expression of AGO, DCL, or RDR homologs [70].

In animals and plants, the presence of conserved core sets of miRNAs helped define the thermodynamic properties for the biogenesis and activity of miRNAs in those systems [10]. However, unlike animals and plants, several recent studies have found little conservation of miRNA genes across fungal taxa [71,72,73,74]. Also, the enzymes involved in RNA silencing are much less conserved in fungi than in higher eukaryotes [75]. For example, neither of the two putative S. sclerotiorum DCLs (SS1G_13747 and SS1G_10369) contained canonical PAZ domains, which could mean that parameters appropriate for identification of pre-miRNAs in S. sclerotiorum may differ markedly from those used for plants. As previously reported [76], algorithms selected for identification of S. sclerotiorum miRNAs and miRNA targets significantly impacted the sequences identified. For example, Zhou et al. [71] described 42 milRNAs in S. sclerotiorum and confirmed that a subset of the milRNAs was differentially expressed during sclerotial development. However, none of the milRNAs reported by Zhou et al. [71] were identified in the present study. This may be due to differences in a priori assumptions of milRNA predictions; 58.5% of the milRNA sequences reported by Zhou et al. [71] were derived from exons, but coding sequences were excluded from our analyses.

Plants and animals possess ancient and recent, i.e., species-specific, miRNAs [77]. Ancient miRNAs have conserved sequences and targets, while species-specific miRNAs have targets that are more diverse. A small subset of the candidate S. sclerotiorum miRNAs was conserved in closely related fungal species as reported for oomycetes in the genus Phytophthora [56]. The most broadly conserved small RNA sequences were derived from predicted S. sclerotiorum tRNAs. Small RNAs derived from tRNAs have been shown to regulate gene expression in humans, oomycetes, and plants, and could be important for RNA silencing and gene regulation in S. sclerotiorum [20,78,79,80]. The detection of phasiRNAs suggests that S. sclerotiorum miRNAs are capable of cleaving long noncoding RNAs that act as templates for the production of phasiRNAs to regulate other genes in trans [13].

Infection of S. sclerotiorum by SsHV2-L was associated with enhanced production of siRNAs from 437 loci that may be analogous to Arabidopsis vasiRNAs that have been associated with widespread silencing of host genes and the establishment of a broad-spectrum antiviral state [26]. In S. sclerotiorum, 19 of the loci that showed enhanced sRNA accumulation in SsHV2-L-infected samples produced phasiRNAs. In plants, the production of phasiRNAs is initiated by miRNA cleavage of an RNA template, which functions as a substrate for siRNA production at discrete 21 or 22-nt intervals through the combined action of DCL4 and RDR6 [81]. In Arabidopsis, phasiRNAs act as negative regulators of nucleotide-binding, leucine-rich repeat plant defense and other genes that can be triggered by a miRNA [82]. While in Drosophila, phasiRNAs are expressed primarily in germline cells and suppress retrotransposons and may serve similar functions in S. sclerotiorum [81].

Recent studies revealed that fungi are capable of exporting sRNA that can play an important role in pathogenesis [83] and the host can export small RNA to silence fungal genes involved in pathogenesis [84]. The wheat-infecting fungus P. striiformis is predicted to use sRNAs to target resistance gene homologs in wheat [14], which is consistent with findings in B. cinerea [85] that sRNAs serve as fungal virulence factors. Some S. sclerotiorum 22-nt sRNAs from the present study were complementary to mRNAs from B. napus, G. max, and H. annuus mRNAs, including mRNAs predicted to encode ferritin-1. In plants, sRNAs of 22-nt can trigger production of secondary siRNAs and rapid silencing of the targeted mRNA [86]. It remains a question whether a subset of the 22-nt RNAs produced by S. sclerotiorum are exported to plant hosts and trigger production of secondary siRNA to downregulate host defense genes, such as ferritin-1 in plants.

We previously demonstrated that SsHV2L infection significantly reduced virulence, changed mycelial growth patterns in 3-day old cultures, and delayed and reduced sclerotia production [31]. In a previous study, we observed that most of the SsHV2L-derived vsiRNAs originated from antisense RNA [87]. However in this study, SsHV2L-derived vsiRNAs were produced from both sense and antisense RNAs, which could be related to the different growth conditions used for the two studies. In summary, we predicted novel classes of sRNAs, showed that infection by a mycovirus induced the expression of sRNAs from both coding and noncoding RNAs, and showed that S. sclerotiorum produced large amounts of tRNA-derived siRNAs and phasiRNAs, the latter of which may be capable of regulating gene expression in trans. The roles of these new classes of fungal sRNAs in gene regulation remain to be determined. Future experiments to dissect fungal RNA silencing pathways by introducing mycoviruses may allow us to assign functions to each putative S. sclerotiorum RNA silencing gene. A better understanding of small RNA processing will allow us to use the obtained knowledge to compromise fungal antiviral defenses to better control fungal pathogens.