FoQDE2-dependent milRNA promotes Fusarium oxysporum f. sp. cubense virulence by silencing a glycosyl hydrolase coding gene expression

MicroRNAs (miRNAs) are small non-coding RNAs that regulate protein-coding gene expression primarily found in plants and animals. Fungi produce microRNA-like RNAs (milRNAs) that are structurally similar to miRNAs and functionally important in various biological processes. The fungus Fusarium oxysporum f. sp. cubense (Foc) is the causal agent of Banana Fusarium vascular wilt that threatens global banana production. It remains uncharacterized about the biosynthesis and functions of milRNAs in Foc. In this study, we investigated the biological function of milRNAs contributing to Foc pathogenesis. Within 24 hours post infecting the host, the Argonaute coding gene FoQDE2, and two Dicer coding genes FoDCL1 and FoDCL2, all of which are involved in milRNA biosynthesis, were significantly induced. FoQDE2 deletion mutant exhibited decreased virulence, suggesting the involvement of milRNA biosynthesis in the Foc pathogenesis. By small RNA sequencing, we identified 364 small RNA-producing loci in the Foc genome, 25 of which were significantly down-regulated in the FoQDE2 deletion mutant, from which milR-87 was verified as a FoQDE2-depedent milRNA based on qRT-PCR and Northern blot analysis. Compared to the wild-type, the deletion mutant of milR-87 was significantly reduced in virulence, while overexpression of milR-87 enhanced disease severity, confirming that milR-87 is crucial for Foc virulence in the infection process. We furthermore identified FOIG_15013 (a glycosyl hydrolase-coding gene) as the direct target of milR-87 based on the expression of FOIG_15013-GFP fusion protein. The FOIG_15013 deletion mutant displayed similar phenotypes as the overexpression of milR-87, with a dramatic increase in the growth, conidiation and virulence. Transient expression of FOIG_15013 in Nicotiana benthamiana leaves activates the host defense responses. Collectively, this study documents the involvement of milRNAs in the manifestation of the devastating fungal disease in banana, and demonstrates the importance of milRNAs in the pathogenesis and other biological processes. Further analyses of the biosynthesis and expression regulation of fungal milRNAs may offer a novel strategy to combat devastating fungal diseases.


Introduction
Banana Fusarium wilt, also known as Panama disease, is caused by the fungal pathogen Fusarium oxysporum f. sp. cubense (Foc). And it poses a serious threat to the banana industry worldwide [1,2]. Four physiological races of Foc have been identified, of which tropical race 4 (TR4) has the most devastating effect on banana production [2,3]. TR4 originated in Southeast Asia and is the main cause of banana Fusarium wilt in China [3][4][5][6]. TR4 has spread rapidly around the world, and has been reported in Mozambique, Australia, Pakistan, and even in countries along the Mediterranean coast, in countries such as Lebanon, Oman, and Jordan [1,2,7]. However, fungicides, flood fallowing, and organic amendments have rarely provided long-term control in any banana planting area. The only effective method for controlling the dissemination and subsequent infections by Foc in banana is by the quarantine or exclusion of infected properties or by planting non-host crops or cultivars [8]. Lack of effective methods to control banana fusarium wilt seriously imperils global banana production. Improved strategies to control this devastating disease are urgently needed.
MicroRNAs (miRNAs), a type of small non-coding single-stranded RNAs, play crucial roles in diverse biological processes [9,10]. Through base pairing with target messenger RNAs (mRNAs), miRNAs degrade the target mRNA or inhibit its translation and thereby regulate gene expression at the post transcription level [11]. Small RNAs (sRNAs) have been reported in various fungi, including Neurospora crassa, Magnaporthe oryzae, Botrytis cinerea, and Sclerotinia sclerotiorum [12][13][14][15]. Because some of these sRNAs are structurally similar to miRNAs from plants and animals, they are called microRNA-like RNAs (milRNAs) [12,14]. In B. cinerea, milRNAs have been identified as virulent effectors that suppress host immunity and facilitate fungal infection [14,16]. In Verticillium dahliae, a novel milRNA, VdmilR1, was reported to play crucial roles in pathogenicity [17]. Recently, a milRNA (Vm-milR37) expressed exclusively in the mycelium, was verified to contribute to pathogenicity in Valsa mali [18]. On the other hand, plants also export miRNAs or siRNAs to inhibit gene expression in fungal pathogens and confer efficient crop protection from pathogen infection [19][20][21]. Thus, the trans-kingdom sRNAs play key roles in host-pathogen interactions [22].
Four types of milRNAs generated from four different biosynthesis pathways, namely milR-1 to -4, have been reported in the model fungus N. crassa [12]. Different combinations of factors, including Dicers, QDE2 (Quelling Deficient 2), the exonuclease QIP (QDE2 interacting protein), and an RNAse III domain-containing protein, MRPL3, are involved in the production of milRNAs [12,23,24]. The reported milRNA biosynthesis pathways in N. crassa appear more complex and diverse than that in plants and animals [25].
Argonaute (AGO) proteins are evolutionarily conserved in all domains of life and play a key role in the RNA interference (RNAi) pathway [26]. As in plants and animals, AGO proteins are the core components of the RNA-induced silencing complex (RISC) and contribute to gene silencing by using the loaded sRNA as a specificity determinant in fungi [15]. QDE2, one of the two identified AGO proteins in N. crassa, functions as a slicer and is required for the biogenesis of some sRNAs such as milRNAs and PIWI-interacting RNAs [12,24]. Suppressor of meiotic silencing 2 (Sms-2), another reported AGO protein in the N. crassa genome, is thought to function by binding to sRNAs originating from the unpaired DNA region and is required for meiotic silencing of unpaired DNA [25,27]. In the model fungus N. crassa, milR-1 type of milRNAs are the most abundant milRNAs and the maturation of the milRNAs requires the AGO protein QDE-2 [12,24]. However, whether this type of milRNAs exist and what function they have in Foc remain unclear.
In this study, we identified an AGO and two Dicer proteins in Foc, and examined their functions in milRNAs biosynthesis and in fungal pathogenesis by sRNA sequencing and reverse genetics. We identified a FoQDE2-dependent milRNA (milR-87), which contributes to invasive growth during the initial stage of Foc infection and thus affects Foc pathogenicity, likely by targeting the gene FOIG_15013. FOIG_15013 encodes a glycosyl hydrolase, which appears as a negative regulator of Foc conidiation and pathogenicity. Overall, our findings uncover the novel function of milRNA in Foc pathogenicity.

Identification of AGO protein FoQDE2 in Foc
By performing an orthologous protein BLAST (Basic Local Alignment Search Tool, https:// blast.ncbi.nlm.nih.gov/Blast.cgi) search, we identified two AGO proteins in Foc, encoded by FoQDE2 (FOIG_01986) and FoAGO2 (FOIG_01246) respectively. Conserved domain prediction showed that the proteins have four domains: a variable N-terminal domain (ArgoN), a linker 1 domain (ArgoL1), a PAZ domain, and a PIWI domain (Fig 1A).
The amino acid sequences of AGO orthologs in other fungal species were retrieved from GenBank for phylogenetic analysis. Two to five AGO-like proteins were identified in different strains of F. oxysporum, but only two in Foc. Phylogenetic analysis indicated that AGO proteins from filamentous fungi could be divided into three subgroups (Fig 1B). In the first group, FoQDE2 and its orthologs from F. oxysporum and F. graminearum were clustered with QDE2 from N. crassa, AGO-like protein MoAGO3 from M. oryzae [15], and Agl1 from the chestnut blight fungus Cryphonectria parasitica [28]. The second group includes AGO-like proteins from F. oxysporum and Agl2 and Agl3 from C. parasitica. The third group is composed of FoAGO2 and its orthologs from F. oxysporum and F. graminearum, Agl4 from C. parasitica, MoAGO1 from M. oryzae, and SMS2 from N. crassa [25]. The phylogenetic analysis showed that the AGO proteins in Foc are orthologous to those of the other filamentous fungi. . Bootstrap values are expressed as a percentage of 1000 replicates. In addition to Argonaute proteins from the Fusarium oxysporum f. sp. cubense (Foc) genome showed in bold, proteins used for this analysis include Agl1-Agl4 from Cryphonectria parasitica [28], a chestnut blight fungus; MoAGO1-MoAGO3 from Magnaporthe oryzae, a rice blast fungus [15]; NcQDE2 and NcSMS-2 from the model fungus Neurospora crassa [25]

Characterization of sRNA biosynthesis-related genes in Foc
In addition to the AGO coding genes, two Dicer coding genes (FoDCL1 and FoDCL2) exist in Foc, which function in the biosynthesis of fungal sRNAs [12]. Compared with Foc in pure culture conditions, FoQDE2 and Dicer coding genes were significantly up-regulated at 24 hours post inoculation (hpi) on host banana plants. In contrast, transcription of FoAGO2 was nearly undetectable (Fig 2A). The transcriptional induction of these sRNA biosynthesis genes during host infection stage suggests that the sRNAs synthesis mediated by these genes may play an active role during Foc pathogenesis.
To further investigate the function of the three up-regulated genes (FoQDE2, FoDCL1, and FoDCL2), we generated these genes deletion mutants in the wild-type (WT) strain XJZ2 via homologous recombination (S1A Fig). Four independent FoQDE2-deletion mutants were verified by PCR and Southern blot analysis (S1B and S1C Fig). Four FoDCL1-deletion mutants and FoDCL2-deletion mutants were respectively obtained and verified by Southern blot analysis (S1C Fig). Four complemented strains were also generated by expressing the  Table).
Relative transcript levels of FoQDE2, FoDCL1, and FoDCL2 were examined by quantitative real-time PCR (qRT-PCR) analysis in the WT strain and the mutants. No FoQDE2, FoDCL1, or FoDCL2 transcript was detected in the corresponding deletion mutants (S2A Fig), confirming that the genes had been successfully deleted in the respective mutant. Transcript level of FoQDE2 was detected in the complemented strains, to a comparable level of that in the WT (S2B Fig), indicating that FoQDE2 was restored in the complemented strains (cΔFoQDE2-1/-2).

Characterization of FoQDE2 during growth, conidial production, stress tolerance and virulence to banana seedlings in Foc
We next assessed the biological function of FoQDE2 by examining growth and conidiation phenotypes in the ΔFoQDE2 mutants. The results showed that, the colonial morphology of the ΔFoQDE2 was strikingly different from that of the WT strain XJZ2. The aerial hyphae of ΔFoQDE2 were much less abundant than those of the WT, and the mycelia grew close to the surface of the PDA media. Such morphological change of the colonies was restored to different degrees in the two complemented FoQDE2 strains (cΔFoQDE2-1 and cΔFoQDE2-2; Fig 2C). On the other hand, the colony morphology was not much changed in the ΔFoDCL1 or ΔFoDCL2 mutants (Fig 2C).
We assessed mycelial growth rate in the WT and mutant strains. From the fourth day of culture, the growth of ΔFoQDE2 was significantly slower than that of the WT at the level of α = 0.01. The growth of the complemented strains was not fully recovered. The ΔFoDCL1 and ΔFoDCL2 mutants displayed no difference in mycelial growth compared with the WT (S2 Table). Thus, loss of FoQDE2 led to slow mycelial growth on PDA medium. Furthermore, microconidia production was significantly reduced in the ΔFoQDE2 mutant, compared to the WT, when cultured on PDA medium for 7 days. Such conidiation defect could be partially restored in the FoQDE2 complemented strains (Fig 2B). Compared with the WT, ΔFoDCL1 also produced less microconidia, whereas ΔFoDCL2 showed no significant difference in conidial production ( Fig 2B).
To elucidate the causes of the decreased conidia production in ΔFoQDE2, six reported conidiation-related genes [29][30][31][32] were selected for their transcript levels assessment by qRT-PCR. Compared with the WT, the transcript levels of four selected conidiation-related genes, StuA,

PLOS PATHOGENS
FoQDE2-dependent milRNA promotes Fusarium oxysporum f. sp. cubense virulence Overall, our results showed that deletion of FoQDE2 leads to morphological changes, including changes in colonial morphology and reduction of mycelial growth, as well as reduction in conidial production of Foc.
We next tested the sensitivities of the WT and the mutant strains to oxidative stress. The ΔFoQDE2 mutant was hypersensitive to 3 mM H 2 O 2 when cultured on minimal medium (MM), which could be restored in the complemented strains (Figs 2C and S2D). However, the ΔFoDCL1 and ΔFoDCL2 mutants showed no difference in sensitivity to oxidative stress as compared to the WT when cultured on MM supplemented with 2 or 3 mM H 2 O 2 (Figs 2C and S2D). Taken together, FoQDE2 is required for tolerance to oxidative stress.
Mycelial infection capacity was tested on the tomato fruit surface and the incidence of tomato fruit infection was recorded at five days post inoculation (dpi). The ΔFoQDE2 mutant did not penetrate the epidermis or cause fruit tissue necrosis on the surface of tomato as the WT and complemented strains did ( Fig 3A). In contrast, the ΔFoDCL1 and ΔFoDCL2 mutants caused tissue necrosis similar to that of WT ( Fig 3A). Thus, FoQDE2 is required for successful invasive growth.
Pathogenicity tests on banana (Cavendish) seedlings showed that the ΔFoQDE2 mutant was unable to produce obvious vascular discoloration in the corm of banana seedling (Fig 3B and  3C) and was significantly reduced in virulence on banana compared to the WT (p<0.01, Fig  3D and 3E). The disease index of the complemented strain (cΔFoQDE2-1) was increased but did not reach the level of the WT (Fig 3E), which exhibited brown discoloration in the corm of banana seedlings. On the other hand, compared to the WT, pathogenicity of the ΔFoDCL1 mutant was not changed, while pathogenicity of the ΔFoDCL2 mutant was decreased significantly ( Fig 3D and 3E). These observations indicate that FoQDE2 is involved in the virulence of Foc.

Identification of FoQDE2-dependent milRNA in Foc
To investigate the function of FoQDE2 in sRNA production, we sequenced sRNAs from the WT strain and ΔFoQDE2 mutant. The sequencing data were deposited in the NCBI sequence read archive (SRA) under accession number PRJNA562097. Origin analysis of the sRNAs showed that reads mapped to tRNA in ΔFoQDE2 was significantly reduced by approximately 40% than in the WT, whereas significantly more reads mapped to the UTRs, intronic, and intragenic regions in ΔFoQDE2 than in the WT (Fig 4A). sRNA length distribution analysis showed that most reads with length of between 16-24 nt were in both WT and ΔFoQDE2 mutant, while reads of 21 and 22 nt were reduced in the ΔFoQDE2 mutant compared to the WT ( Fig 4B). Reads starting with A or U were more abundant than those starting with C or G ( Fig 4C). To identify sRNA-producing loci, reads with counts of 10 or higher that mapped to the UTR and intronic and intergenic regions were analyzed. Overlapping and adjacent reads were grouped and only those with a consensus length of less than 300 nt were considered to be small RNA-producing loci. Using this method, 364 loci were captured in the WT and ΔFoQDE2 mutant. Among them, 25 loci were significantly decreased and 13 loci were significantly increased in ΔFoQDE2 compared with the WT. (Fig 4D). A total of 12 sRNAs coming from these loci were selected for expression verification through reverse transcription-PCR.
And the results of electrophoresis showed that the expression of selected sRNAs was consistent with the results of sequencing when using small nuclear RNA U4 of Foc as the internal control (S3A and S3B Fig).
To identify FoQDE2-dependent milRNA, sRNA read counts obtained from the WT and ΔFoQDE2 were normalized and compared. To eliminate interference from sRNAs produced by mRNA degradation, the reads mapped to ORFs were discarded. A total of 25 sRNA-

PLOS PATHOGENS
FoQDE2-dependent milRNA promotes Fusarium oxysporum f. sp. cubense virulence producing loci were significantly decreased in ΔFoQDE2 compared with the WT. Sequences of the loci were extracted from the Foc II5 genome and uploaded to the RNAfold web server (http://rna.tbi.univie.ac.at//cgi-bin/RNAWebSuite/RNAfold.cgi) for RNA advanced structure prediction. The loci with a stem-loop RNA structures were validated as precursors of milR-NAs. Among them, four milRNAs were predicted to form stem-loop structures. Three of them

PLOS PATHOGENS
FoQDE2-dependent milRNA promotes Fusarium oxysporum f. sp. cubense virulence (milRNA loci were marked as INFOIG_12093, DSFOIG_07638 and DSFOIG_15767) were confirmed with down-regulated expression in the ΔFoQDE2 mutant by reverse transcription-PCR (S3A Fig). We failed to detect expression of the fourth one (UTRFOIG_14681) for its rich in Adenine and Uracil nucleotides. In this study, we focused on one milRNA (milR-87 located at the downstream of FOIG_15767, DSFOIG_15767) that formed a stem-loop precursor ( Fig  5A). The mature milR-87 was predicted to be of 23 nt (Fig 5A). T-clone and sequencing of milR-87 confirmed the results. We further verified expression level of milR-87 using qRT-PCR ( Fig 5B) and sRNA Northern blot (S4E Fig). Compared to the WT, the production of milR-87 was significantly reduced in the ΔFoQDE2 and ΔFoDCL2 mutants (Figs 5B and S4E). Based on our sRNA analysis, milR-87 is dependent on FoQDE2 and FoDCL2 in Foc.

MilRNA (milR-87) regulates tolerance to hydrogen peroxide and is involved in the virulence of Foc
To determine the function of the milR-87, we generated deletion or overexpression mutants of The ΔmilR-87 mutant was hypersensitive to the oxidative stress generated by 3 mM H 2 O 2 , a phenotype similar to that of the ΔFoQDE2 mutant ( Fig 5D). In contrast, the growth of milR-87 overexpress transformant (OEmilR-87-1) was not affected by oxidative stress (Fig 5D). The result indicates milR-87 positively regulates tolerance to oxidative stress in Foc.
Pathogenicity tests showed that the ΔmilR-87 mutant was compromised in penetrating the cellophane membrane and caused slighter necrosis on the surface of tomato fruits than the WT. While overexpressed transformant showed the opposite (Fig 5E and 5F). Infection assay using banana seedlings showed that ΔmilR-87 significantly reduced in virulence compared to the WT, while the virulence of overexpressed transformant was enhanced, as it caused obvious internal disease symptoms of brown discoloration (Fig 5G and 5H). Furthermore, qRT-PCR analysis showed that milR-87 was significantly up-regulated during the early stages of infection (24,48, and 96 hpi) compared with levels at 0 hpi conditions (Fig 5C). These results suggest that milR-87 contributes to the infective growth and virulence of Foc.
To confirm the function of milRNA (milR-87) in Foc, synthetic double-strand siRNAs (siRNA-1 and siRNA-2) and single-strand antisense small RNA (inhibitor) were designed (S1 Table) to silence the expression of milR-87. We observed that when the siRNAs or inhibitor was transfected into Foc protoplast inoculated on the banana leaf, the size of the lesion decreased significantly compared to those caused by the Foc protoplasts with water as control (CK) (Fig 6A and 6B). An unrelated single-strand sRNA served as the negative control (NC) did not interfere with lesion formation when added into the Foc protoplasts inoculated on the banana leaf (Fig 6A and 6B). We noticed that the ΔmilR-87 mutant produced the smallest lesion on the banana leaf (Fig 6A and 6B). The qRT-PCR results confirmed that expression of milR-87 could be silenced by the siRNA-1, siRNA-2 and inhibitor to different extents ( Fig 6C). Overall, these results confirmed that milR-87 is required for full virulence of Foc.

MilRNA (milR-87) regulates the virulence by targeting a glycosyl hydrolase gene (FOIG_15013) during the pathogenesis of Foc
To investigate regulatory mechanism of milR-87 in Foc pathogenesis, potential target genes were computationally predicted, and verified by qRT-PCR using total RNAs from the mutants (ΔmilR-87 and ΔFoQDE2) and WT. Totally nine genes were significantly up-regulated in the

PLOS PATHOGENS
FoQDE2-dependent milRNA promotes Fusarium oxysporum f. sp. cubense virulence two mutants compared to the WT (S5 Fig). One of these genes, FOIG_15013, encodes a glycosyl hydrolase, consisting of a signal peptide and a conserved domain of GH-79C (Fig 7A), was selected for further investigation. The predicted targeted site of milR-87 was located in the ORF region of FOIG_15013, we then generated a fusion protein by directly ligating the GFP coding sequence to the C'-terminus of FOIG_15013 (Fig 7A). A point-mutated FOIG_15013 that could not be paired with milR-87, FOIG_15013m, was also fused with GFP (Fig 7A). The constructs were transformed into the ΔmilR-87 mutant or the WT strain respectively, to verify the possible of milR-87 regulating the FOIG_15013 expression by evaluating GFP signal brightness. In the ΔmilR-87 mutant, FOIG_15013-GFP expression was not inhibited, as

PLOS PATHOGENS
FoQDE2-dependent milRNA promotes Fusarium oxysporum f. sp. cubense virulence indicated by the strong GFP signal (Fig 7B and 7C). In contrast, the GFP signal was hardly detected in the WT strain, likely due to milR-87 mediated suppression of FOIG_15013-GFP expression (Fig 7B and 7C). The FOIG_15013m-GFP expression was not affected in the WT background, as milR-87 was unable to pair with the mutation site. Correspondingly, transcript level of FOIG_15013 and its mutated version in different strains was consistent with that of GFP signal detection (Fig 7D), further confirming that FOIG_15013 is targeted by milR-87 in Foc.
In order to clarify whether milR-87 regulates Foc pathogenesis through modulating FOIG_15013 expression, we generated the ΔFOIG_15013 mutants and corresponding complemented transformants (S6 Fig) and characterized them in growth and pathogenicity. The results showed that, the ΔFOIG_15013 mutants grew faster and produced significantly more conidia than the WT strain XJZ2 when cultured on PDA plate. While the gene complemented transformants restored the phenotype (Fig 8A and 8B). Pathogenicity tests showed that the ΔFOIG_15013 mutants were more virulent than the WT, while the complemented transformants exhibited a pathogenicity which is similar to that of the WT (Fig 8C-8F). These results indicate that milR-87 may promote Foc pathogenicity via suppression of FOIG_15013 expression, which is a negative regulator of Foc mycelial growth and virulence.
Hypersensitive response (HR), a form of programmed cell death at the site of pathogen infection, is an important symbol of plant immune activation [33]. To test the effect of FOIG_15013 on host plant immunity, we transiently expressed the milR-87 and its target gene FOIG_15013 in the leaves of Nicotiana benthamiana following the established method [18,33]. Compared to the GFP control, transient expression of FOIG_15013 caused obvious leaf yellowing in N. benthamiana at 48 hours post-agroinfiltration (hpa), while positive control of INF1 caused typical leaf necrosis at the same condition (Fig 8G). And expression of FOIG_15013 could induce significant up-regulation of some pathogenesis-related marker genes including the PR-1, PAL4, LOX, and Osmotin coding gene [34,35], but not induce upregulation of HR-related gene the HIR1 [36] and ethylene response gene the ERF1 [37]. In contrast, co-expression of milR-87 with FOIG_15013 significantly down-regulated these marker genes expression (Fig 8H). The results indicate the FOIG_15013 activates general defense responses in N. benthamiana, while milR-87 suppresses these defense responses likely by silencing the expression of FOIG_15013.

Discussion
Most fungal genomes encode more than one AGO proteins. The N. crassa genome has two, M. oryzae three, C. parasitica four, and members within the F. oxysporum species complex produce two to five AGO proteins [12,15,28]. Eucaryotic AGO proteins are characterized by a bilobal architecture, one contains the N-terminal and PAZ domains, and the other contains the MID and PIWI domains [38]. The Foc genome encodes two AGO proteins, both containing conserved PAZ and PIWI domains. Phylogenetic analysis (Fig 1B) showed the FoQDE2 proteins derived from different formae specialis were closely clustered in one group, suggesting that FoQDE2 proteins are well conserved and the FoQDE2 is orthologous to QDE2 proteins of N. crassa, F. graminearum, and M. oryzae.
In F. oxysporum, several genes related to conidiation were reported [29][30][31][32]. Here we selected six of them to verify the function of FoQDE2 in conidial production. qRT-PCR analysis showed that, except for brlA and Htf1, the transcription of the other four conserved genes (abaA, wetA, FoNIIA, and StuA) related to conidiation in Fusarium or other fungi was downregulated in ΔFoQDE2. The transcript levels of only two of these genes (FoNIIA and StuA) were restored to the WT level in the complementation transformant. Previous studies have

PLOS PATHOGENS
FoQDE2-dependent milRNA promotes Fusarium oxysporum f. sp. cubense virulence shown that a nitrite reductase coding gene, FoNIIA, is significantly up-regulated during conidiation compared with during vegetative growth in F. oxysporum, likely under regulation of a transcription factor Ren1 [29]. StuA is involved in the formation of macroconidia in F. oxysporum and is also related to spore development and pathogenicity in F. graminearum [30,39]. In our phenotypic analysis, conidiation of ΔFoQDE2 was markedly lower than that of the WT, suggesting that FoQDE2 could affect conidial production through down-regulating FoNIIA and StuA expression. In Aspergillus and other fungi, three tandem genes (brlA, abaA, and wetA) are considered to be central regulators of conidial production [31,32]. In this study, the transcript levels of two genes (abaA and wetA) were increased in the complementation transformant, but were not restored to the WT level, suggesting that the decreased conidiation of ΔFoQDE2 was not directly related to this central regulatory pathway. It is unknown whether the reason for the reduced conidiation of the ΔFoQDE2 mutant is related to the FoQDE2-dependent milRNAs.
Previous research showed that the AGO protein FgAgo1 and Dicer protein FgDicer2 function in RNAi-mediated gene silencing in F. graminearum [40]. Deletion of either of the two genes did not affect the growth and pathogenicity of F. graminearum [40]. Our results demonstrate that, compared with WT, the ΔFoDCL1 and ΔFoDCL2 mutants have no differences in the tested phenotypic traits, including mycelial growth and sensitivity to H 2 O 2 . The results are consistent with those reported in F. graminearum [41,42], but not the same as in Va. mali [43]. In our study, the phenotype of the ΔFoQDE2 is completely different from that of F. graminearum. The differences in pathogenicity of Foc and F. graminearum may be caused by differences in infection sites. Foc infects the host through roots, whereas F. graminearum damages the host panicles.
Currently five types of milRNA biosynthetic pathways have been reported in fungi, four reported in N. crassa and one reported in Ve. dahliae [12,17]. Among them, milR-1 and milR-2 produced by the two biosynthetic pathways, respectively, are dependent on the AGO protein QDE2 in N. crassa [12]. Maturation of milR-1 milRNAs, the most abundant milRNAs, requires the AGO protein QDE-2, Dicer, and QIP [24]. In this study, transcript level of milR-87 was significantly lower in ΔFoQDE2 and ΔFoDCL2 mutant, compared with that of the WT. The result indicates milR-87 is a FoQDE2-dependent milRNA and belongs to the milR-1 type of milRNA.
We identified a target gene, FOIG_15013, of milR-87. FOIG_15013 encodes a glycosyl hydrolase, shared high sequence similarity with 4-O-α-L-rhamnosyl-β-D-glucuronidase (FoBGlcA), a novel enzyme of F. oxysproum [44] and β-D-glucuronidase (AcGlcA79A) of Acidobacterium capsulatum [45]. These hydrolytic enzymes are secreted during plant-pathogen interactions and could hydrolyze D-glucuronic acid (GlcA) or 4-O-Me-GlcA residues at the Arabinogalactan proteins (AGPs) side-chain terminals [46,47]. It has been proposed that Arbinogalactans (AGs) are structural components of the plant cell wall and they are associated with proteins forming AGPs. Degradation products of AGPs generated by hydrolytic enzymes of pathogens could function as damage-associated molecular patterns (DAMPs) eliciting the plant defense response [48]. Therefore we speculate that product of FOIG_15013 may secrete out and act as an avirulence effector, and likely activates the host defense responses during the infection of Foc. Pathogenesis-related proteins (PRs) of plant are usually activated by different biotic and abiotic stresses, including pathogen infection, to confer plant resistance [49,50]. Examples of PRs include the PR-1 and PAL4, involved in the salicylic acid (SA)-related defense pathways in N. benthamiana [50,51], the Osmotin belonging to PR-5 family and involved in the plant defense against various pathogens [52], and 9-Lipoxygenase (9-LOX) involved in activating local defense against pathogen [35,53]. In this study, heterologous transient expression of FOIG_15013 in N. benthamiana induced significant up-regulation of these aforementioned resistant marker genes. In contrast, co-expression of milR-87 suppressed the resistant

PLOS PATHOGENS
FoQDE2-dependent milRNA promotes Fusarium oxysporum f. sp. cubense virulence marker genes' expression. These results support that the FOIG_15013 encoding glycosyl hydrolase of GH-79C family may act as an avirulence effector, thus negatively regulates Foc pathogencity, and the milR-87 could suppress the expression of FOIG_15013 to facilitate infection to the host plant. Future investigation is needed to further elucidate the hydrolase activity of FOIG_15013 and interaction between the enzyme and plant resistant marker genes expression.
Some milRNAs, found in phytopathogenic fungi Ve. dahliae and Va. mali, regulate the fungal pathogenicity by targeting their own virulent genes. For example, the VdmilR1, which is independent of Dicer and AGO proteins, regulates fungal virulence at the later stage of inoculation by suppressing a virulence gene (VdHy1) expression through increasing histone H3K9 methylation [17]. A Vm-milR37, exclusively expressed in the mycelium, is involved in pathogenicity by targeting a glutathione peroxidase coding gene VmGP, which contributes to the oxidative stress response [18]. In the study, a milRNA (milR-87), firstly identified in Foc through sRNA sequencing, contributes to the fungal virulence during the early infection stages (24-96 hpi) by targeting its own glycosyl hydrolase coding gene. So it is different from the reported milRNA found in Ve. dahliae and Va. mali.
Previous studies showed that aggressive fungal pathogens such as Botrytis and Verticillium can take up external sRNAs and double-stranded RNAs (dsRNAs) [19,54]. Applying sRNAs or dsRNAs that target Botrytis milRNA synthetic genes BcDCL1 and BcDCL2 on the surface of fruits, vegetables, and flowers significantly inhibits gray mold disease [19]. In our study, the virulence of Foc to banana leaves was significantly reduced by using exogenous siRNAs, a kind of dsRNAs, to suppress and interfere with the synthesis of milR-87, resulting in reduced disease lesion formation mimicking loss of miR87 function. This result suggest that milR-87 could be used as an effective anti-fungal target for developing the new and environment-friendly fungicides against banana Fusarium wilt.
This study provided the undisputable evidence to demonstrate the contribution of milR-87 in promoting the virulence of Foc during the early infection stage. Furthermore, a novel glycosyl hydrolase acting as an avirulence effector was identified as the target of milR-87 to activate general host defense response. And to our knowledge, it is the first report of avirulence effector targeted by fungal milRNA. The results indicate that the milRNA could silence effector coding gene expression to evade the activation of host defense response and ingeniously regulate pathogenicity of F. oxysporum f. sp. cubense.

Fungal strains and culture conditions
The F. oxysporum f. sp. cubense tropical race 4 strain XJZ2, isolated from Guangdong Province in China [55], was used as the WT for fungal transformation, gene knockout, and milRNA overexpression experiments. All fungal strains were cultured on potato dextrose agar (PDA) for conidiation and mycelial growth. Conidiation was induced as described in our previous study [56]. After a 5-day culture on PDA plates, the colony morphology of the tested strains was recorded. Each experiment was repeated three times independently.

Oxidative stress test
To test the sensitivity of strains to hydrogen peroxide, a conidial suspension of each strain was spotted on MM or MM supplemented with H 2 O 2 (2 to 3 mM) and cultured at 28˚C. The colony diameter was measured after 4 days' incubation.

Generation of deletion mutants and gene complementation transformants
Deletion mutants of three RNAi-related genes (FoQDE2, FoDCL1, and FoDCL2) in Foc were created using a conventional target gene replacement method through homologous recombination. For FoQDE2 complementation, a 1,713-bp fragment including the promoter region and ORF of FoQDE2 was cloned into the backbone vector pEX-Zeocin. The ΔFoQDE2-6 mutant was transformed with the complementation vector by PEG-mediated transformation as described previously [56]. For FOIG_15013 gene complementation, the fragment with promoter and ORF was firstly cloned into the vector pKNT-G418 (donated by Dr. YZ Yun from Fujian Agriculture and Forestry University), then was transformed into the ΔFOIG_15013-19 mutant.

Nucleic acid manipulation, small RNA detection and qRT-PCR analysis
Fungal genomic DNA was extracted according to a previously described method [57]. Deletion mutants were verified by PCR and Southern blot analysis according to a previous method [58]. Primers used in the study are shown in S1 Table. Total RNA was extracted using Trizol reagent (Invitrogen, USA) according to the manufacturer's protocol. To detect target gene expression, qRT-PCR was carried out as described previously [57]. At least three independent experiments were performed by using three biological replicates.
For detection of small RNAs, a total of 40-60 μg total RNA was separated in a 15% ureapolyacrylamide gel, and transferred to Amersham Hybond-N + membrane (GE Healthcare, USA). The probes were labelled with biotin using Biotin 3´End DNA Labeling Kit (Thermo scientific, USA). Hybridization signals were detected by Chemiluminescent Nucleic Acid Detection Module (Thermo scientific, USA). Signal intensity was quantified using Image Lab 6.1.0 software (BIO-RAD, USA).
An All-in-One miRNA qRT-PCR Detection Kit (GeneCopoeia, Rockville, MD) was used for milRNA expression analysis [13]. Briefly, a PolyA tail was added to total RNAs using PolyA polymerase (NEB, USA). The RNA was then reverse-transcribed with an oligo-dT adaptor. MilRNA expression was detected using qRT-PCR, as described above, except that the reference gene was replaced by snRNA U4 in Foc. At least three biological replicates for each sample were performed.
For small RNA expression verification, reverse transcription-PCR were used in the study [42]. Briefly, cDNA was reverse-transcribed from the RNAs added with a PolyA tail with an oligo-dT adaptor. Reverse transcription-PCR was set with 20-26 cycle to avoiding nonspecific amplification.

Small RNA library construction and sequencing
Total RNAs of the WT strain XJZ2 and FoQDE2 deletion mutant ΔFoQDE2 were extracted as described above. RNA integrity was assessed and quantified with a Bioanalyzer 2100 (Agilent, USA) and only qualified RNA samples (RNA integrity number, RIN: 7-10) were used for sRNA library construction. Three different samples of every strain were used for sRNA library preparation and sRNA sequencing. Small RNA libraries were prepared using NEB Next Small RNA Library Prep Set for Illumina (NEB, USA) according to the manufacturer's protocol. Briefly, total RNAs were reverse-transcribed and indexed. Then, cDNAs were separated on a 6% polyacrylamide gel and cDNAs ranging from between 140 and 160 nt, corresponding to 20-40 nt sRNAs, were cut and recycled. The size and concentration of the sRNA sequencing libraries were assessed again using Bioanalyzer 2100. Three sRNA libraries from the WT fungus and two from the ΔFoQDE2 mutant were qualified for sRNA sequencing. Sequencing was performed on an Illumina MiSeq platform at the University of Massachusetts, Amherst using MiSeq Reagent Kit v2 for single-end reads at a length of 50 bp.

Small RNA data analysis and milRNA prediction
Raw sequencing data were trimmed and analyzed using CLC Genomics workbench software v10.1 (CLC bio) and only reads with a read count number larger than 10 and length between 16 and 40 nt were kept for further analysis. The Foc II5 (TR4 strain) genome released by the Broad Institute was used as a reference sequence. To examine the origin of the sRNA, reads mapped to the FocII5 genome were collected and sequentially mapped to the rRNA, tRNA, snRNA, snoRNA, repeat, exonic, intronic, and 3 0 and 5 0 UTRs regions of coding genes, and intergenic regions (the region 1000 nt upstream and 1000 nt downstream of genes). To reduce false positive sRNAs induced by mRNA degradation, only reads mapped to the UTR, intron, and intergenic regions were used for sRNA length distribution and sRNA locus identification.
Overlapping and adjacent reads were grouped and only sequences with a consensus length of less than 300 nt were selected and uploaded to the RNAfold web server (http://rna.tbi.univie. ac.at//cgi-bin/RNAWebSuite/RNAfold.cgi) for RNA secondary structure prediction.
To identify FoQDE2-dependent sRNAs and their loci, reads from the WT and ΔFoQDE2 were compared to each other. Read counts were normalized by calculating their transcripts per million (TPM) value. After normalization, log2 ratios were calculated using the TPM value from two ΔFoQDE2 mutants and three WT XJZ2 strains. Only sRNAs with log2 ratios of less than -1 or greater than 1 and a p value of less than 0.05 were considered as significantly expressed. The sRNAs, which were significantly down-regulated in ΔFoQDE2 and their precursors with a typical secondary stem and loop structure were considered as FoQDE2-dependent milRNAs. For all predicted and differentially expressed milRNAs, transcript levels were verified by reverse transcription-PCR analysis.

Generation of the milRNA deletion mutants and overexpression transformants
The milRNA deletion mutants were created by replacing the precursor (about 300 bp) of milRNA through homologous recombination as functional gene deletion we previously described [56]. The vector pSilent-1 was donated by Dr. Hitoshi Nakayashiki for milRNA overexpression. The milRNA precursor was amplified from the Foc genome and then cloned into the pSilent-1 vector through digestion with HindIII and KpnI. The inserted fragment was verified by sequencing and the correct vector was used for transformation of the WT strain XJZ2 of Foc. Transformants were identified by PCR with primers given in the S1 Table.

Transfection of siRNAs and inhibitor to suppress the production of milR-87 in Foc
Protoplast preparation was performed as described previously [56]. Double-strand siRNAs (siRNA-1 and siRNA-2) and single-strand inhibitor against the precursor sequence of milR-87, as well as a random sRNA sequence (24 nt) serving as a negative control (NC) were designed and synthesized by General Biosystems company. According to the protocol provided by transfection reagent, a total of 100 nM siRNA oligos or inhibitor was transfected into the prepared protoplast of Foc using lipofectamine 2000 (Invitrogen). After two days culture, conidia suspensions of the transfected strains and positive control of the ΔmilR-87 mutant, as well as non-transfected strain of CK were inoculated on the banana leaves. After three days, the lesions on banana leaves were recorded and measured for accessing virulence of different treatments. The total RNAs of transfected strains cultured at 28˚C for 7 days were extracted for milRNA detection.

Fungal invasion assays and pathogenicity test
The ability of mycelia to penetrate cellophane and the invasive growth on tomato fruit surfaces were compared between WT and the mutants according to a previous description [56]. The pathogenicity on banana seedlings was assessed as described previously [5]. A total of 30 banana plantlets were used for each treatment. The severity of internal symptoms was recorded according to the disease grade [5], and disease indexes were calculated for disease severity assessment.

Statistical analysis
In this study, a Student's t-test was used for significance analysis of two samples. While a Duncan's multiple range test was used to analyze the data of multiple samples at the both levels (α = 0.05 and α = 0.01). And different letters indicated significant difference.

MilRNA target gene identification and transient expression in Nicotiana benthamiana
Based on the sequence of milR-87, its target genes in Foc genome were predicted using psRNA-Target online software. In which the genes significantly up-regulated in the mutants ΔmilR-87 and ΔFoQDE2 were screened by qRT-PCR and selected for target gene identification. According to the target site of the milRNA, a GFP marker gene was fused to the 3 'terminal of the candidate target gene. And a point-mutation of the milR-87 targeted site was introduced into the WT. The GFP fluorescence intensity representing the expression of the target gene was accessed by confocal microscopy. Furthermore, the expression of GFP and target genes was quantified by qRT-PCR.
Transient expression vector was constructed according to the method described by [33]. Briefly, precursor of the milR-87 and FOIG_15013 were introduced to vector pBIN:GFP. The recombinant vectors were transformed into Agrobacterium tumefaciens GV3101. For transient expression, transformed A. tumefaciens cultures were injected into N. benthamiana leaves [14]. After 48 h, the injected leaves were photographed and then harvested for detecting mRNA and protein levels of the FOIG_15013, as well as transcript level of resistant marker genes described previously in N. benthamiana [35,36,[50][51][52].
Supporting information S1