Genome-Wide Identification of Dicer-Like, Argonaute, and RNA-Dependent RNA Polymerase Gene Families in Brassica Species and Functional Analyses of Their Arabidopsis Homologs in Resistance to Sclerotinia sclerotiorum

RNA silencing is an important mechanism to regulate gene expression and antiviral defense in plants. Nevertheless, RNA silencing machinery in the important oil crop Brassica napus and function in resistance to the devastating fungal pathogen Sclerotinia sclerotiorum are not well-understood. In this study, gene families of RNA silencing machinery in B. napus were identified and their role in resistance to S. sclerotiorum was revealed. Genome of the allopolyploid species B. napus possessed 8 Dicer-like (DCL), 27 Argonaute (AGO), and 16 RNA-dependent RNA polymerase (RDR) genes, which included almost all copies from its progenitor species B. rapa and B. oleracea and three extra copies of RDR5 genes, indicating that the RDR5 group in B. napus appears to have undergone further expansion through duplication during evolution. Moreover, compared with Arabidopsis, some AGO and RDR genes such as AGO1, AGO4, AGO9, and RDR5 had significantly expanded in these Brassica species. Twenty-one out of 51 DCL, AGO, and RDR genes were predicted to contain calmodulin-binding transcription activators (CAMTA)-binding site (CGCG box). S. sclerotiorum inoculation strongly induced the expression of BnCAMTA3 genes while significantly suppressed that of some CGCG-containing RNA silencing component genes, suggesting that RNA silencing machinery might be targeted by CAMTA3. Furthermore, Arabidopsis mutant analyses demonstrated that dcl4-2, ago9-1, rdr1-1, rdr6-11, and rdr6-15 mutants were more susceptible to S. sclerotiorum, while dcl1-9 was more resistant. Our results reveal the importance of RNA silencing in plant resistance to S. sclerotiorum and imply a new mechanism of CAMTA function as well as RNA silencing regulation.


INTRODUCTION
RNA silencing refers to a variety of mechanisms whereby a small RNA molecule interferes with a given nucleotide sequence. It was first discovered in plants and occurs widely in eukaryotic organisms (Tijsterman et al., 2002;Ullu et al., 2004). In plants, RNA silencing is triggered by double-stranded RNA (dsRNA) that gives rise to small RNAs known as microRNAs (miRNAs) or small-interfering RNAs (siRNAs). Generation and function of these small RNAs depend on three key protein families, Dicerlike proteins (DCLs), Argonautes (AGOs), and RNA-dependent RNA polymerases (RDRs; Baulcombe, 2004). A whole RNA silencing process comprises three stages: initiation, maintenance, and signal amplification. DCLs undergo RNase III-type activities to process complementary double-strand RNAs into small RNAs with 21-26 nucleotides in length (Carmell and Hannon, 2004). These small RNAs are then incorporated into AGO-containing RNA-induced silencing complexes (RISCs) to serve as the sequence specificity in RNA degradation, translational inhibition, or heterochromatin formation (Bologna and Voinnet, 2014). At the signal amplification stage, RDR enzymes are responsible for synthesis of dsRNAs from ssRNA templates to initiate a new round of RNA silencing (Sijen et al., 2001).
DCL, AGO, and RDR are key components of RNA silencing machinery. DCL proteins are key components in small RNA biogenesis. They are characterized by the presence of six domains: DEAD, helicase-C, DUF283, PAZ, RNase III, and dsRNA-binding motif (DSRM; Margis et al., 2006). DCL consists of a small gene family in higher plants, insects, protozoa, and some fungi, whereas only one Dicer protein exists in vertebrates, nematodes, Schizosaccharomyces pombe, and green alga Chlamydomonas reinhardtii (Liu et al., 2009). AGO proteins are core factors of the RISC that guide small RNAs to their targets by sequence complementarity, which results in target mRNA cleavage, translational repression, or chromatin modification (Hannon, 2002;Moazed, 2009). AGOs are large proteins (∼90-100 kDa) consisting of several characteristic functional domains including DUF1785, PAZ, MID, and PIWI (Hutvagner and Simard, 2008). RDRs enhance the potency of RNAi by amplifying the aberrant RNA population. It is defined by the presence of a conserved RNA-dependent RNA polymerase catalytic domain and is required for initiation and amplification of the silencing signal (Schiebel et al., 1998). Multiple copies of AGO and RDR genes are known to exist in both plants and animals. Members of these gene families play different roles in RNA silencing. For example, the Arabidopsis thaliana genome contains four DCL proteins (DCL1-DCL4) that specifically produce different types and sizes of small RNAs (Bologna and Voinnet, 2014). The role of DCL1 is to extract a single small RNA duplex out of a RNA loop called the pri-miRNA. This gives rise to a miRNA, generally 21 nucleotides long, and is typically involved in regulating developmental genes (Parent et al., 2015). DCL2 can produce abundant 22-nt viral siRNAs and shares functional overlap with DCL4 in antivirus defense (Moissiard et al., 2007). DCL3 generates 24-nt repeat-associated siRNAs (ra-siRNAs) and is involved in antiviral defense against DNA viruses (Moissiard and Voinnet, 2006). DCL4 generates 21-nt trans-acting siRNAs (ta-siRNA) and is the primary DCL component of antiviral defense against RNA viruses (Deleris et al., 2006). Likewise, AGO1, the most well-studied plant AGO gene, associates with miRNAs and some siRNAs such as ta-siRNAs to cleave target mRNA and/or inhibit translation (Yu and Wang, 2010). AGO2 protein is involved in antiviral defense by catalyzing viral RNA cleavage in Arabidopsis (Jaubert et al., 2011). AGO10, the closest paralog of AGO1, is functionally redundant with AGO1 in some aspects of development (Lynn et al., 1999) and also functions as a decoy for miR165/166 to prevent the formation of AGO1-miR165/166 complexes and the subsequent repression of HDZIP III gene expression . For RDRs, RDR2 converts ssRNAs generated from repetitive DNAs to precursor dsRNAs of ra-siRNAs , while RDR6 produces the ta-siRNA precursors (Yoshikawa et al., 2005). Gene families encoding these three key components of RNA silencing machinery have been identified only in several plant species such as A. thaliana, Oryza sativa (Kapoor et al., 2008), Zea mays (Qian et al., 2011), Solanum lycopersicum (Bai et al., 2012), Nicotiana benthamiana (Nakasugi et al., 2013), Setaria italica (Yadav et al., 2015), and Vitis vinifera (Zhao et al., 2015). Identification of these families in more plant species will enhance our understanding of RNA silencing.
RNA silencing plays multiple roles in regulating growth and development as well as abiotic and biotic stress responses. In higher plants, RNA silencing functions as an antiviral defense through the action of DCL, AGO, and RDR proteins (Ding and Voinnet, 2007). The importance of RNA silencing in plant viral defense is manifested by the fact that it has elicited counter defense measures from the viral pathogens to overcome it. Plant viruses have evolved various viral RNA silencing repressors (VSR) to counteract this defense mechanism by targeting different RNA silencing pathway components (Csorba et al., 2015). Apart from viral defense, evidence accumulates for RNA silencing to play a role in plant interactions with bacterial pathogens (Voinnet, 2008). The first example is a miRNA from Arabidopsis that contributes to basal defense against Pseudomonas syringae by regulating auxin signaling (Navarro et al., 2006). Similar to viruses, the bacteria has also developed mechanisms to suppress RNA silencing in order to cause disease (Navarro et al., 2008). Recently, through the use of mutants for key components of RNA silencing or functional analyses of miRNAs in plant defense, the potential role of RNA silencing in plant defense against fungal pathogens has been revealed. These fungi include Verticillium dahliae (Ellendorff et al., 2009), Verticillium longisporum (Shen et al., 2014), Magnaporthe oryzae (Li et al., 2014), and Botrytis cinerea (Jin and Wu, 2015).
Brassica napus is an allotetraploid and was formed about 7500 years ago by crossing between B. oleracea and B. rapa, followed by chromosome doubling (Chalhoub et al., 2014). It is one of the most important oil crops, yet few RNAi machinery components have been characterized to date. We have identified the miRNAs involved in the interactions between B. napus and Sclerotinia sclerotiorum, one of the most devastating fungal pathogens in oil and vegetable crops. Furthermore, we find that Arabidopsis ago1 and ago2 mutant plants exhibit enhanced susceptibility to S. sclerotiorum (Cao et al., 2016). Our results provide a clue to the important roles of RNA silencing in the interactions between B. napus and S. sclerotiorum. In this study, taking advantage of the completion of the B. napus genome sequencing (Chalhoub et al., 2014), we performed comprehensive bioinformatics analyses to identify DCL, AGO, and RDR gene families that are the three key components of RNA silencing machinery in B. napus. Furthermore, we employed mutants to probe their functions in resistance to S. sclerotiorum. We revealed the significant difference in RNA silencing machinery composition between B. napus and Arabidopsis, demonstrated the important role of RNA silencing in resistance to S. sclerotiorum and indicated a possible regulating mechanism of RNA silencing.

MATERIALS AND METHODS
Identification of Putative B. napus DCL, AGO, and RDR Genes Protein sequences of Arabidopsis DCLs, AGOs, and RDRs were downloaded from TAIR database (http://www.arabidopsis.org/) and scan for conserved domains were performed using National Center for Biotechnology Information Conserved Domain Database (NCBI-CDD; http://www.ncbi.nlm.nih.gov/Structure/ cdd/wrpsb.cgi). All these protein sequences were used as queries to search their orthologs in B. napus, B. rapa, and B. oleracea genomes using BLASTp program in NCBI database with default settings. All retrieved protein sequences were examined for the presence of conserved domains and redundant sequences were removed. All candidate sequences of B. napus were subsequently verified in the GENOSCOPE database (http://www.genoscope. cns.fr/blat-server/cgi-bin/colza/webBlat). The physico-chemical properties of BnDCL, BnAGO, and BnRDR proteins were then predicted using ExPASy Compute pI/Mw tool (http://au.expasy. org/tools/pi_tool.html; Bjellqvist et al., 1993).

Phylogenetic Analysis and Nomenclature
Multiple alignment of DCL, AGO, and RDR protein sequences from A. thaliana, B. napus, B. rapa, and B. oleracea was performed using MUSCLE program (Edgar, 2004). The phylogenetic trees were then constructed using MEGA 5.0 (Tamura et al., 2011) by Neighbor-Joining (NJ) method following the Jones-Taylor-Thornton (JTT) model. Bootstrap analysis was performed with 1000 replicates to assess statistical support for nodes. The candidate genes were renamed according to the phylogenetic relationship and sequence homologies with corresponding A. thaliana homologs. The detail information of the proteins used for phylogenetic tree construction was listed in Table 1 and  Table S1.

Plant Materials and Inoculation Analyses
The Arabidopsis dcl, ago, and rdr mutants were provided by Prof. Shou-Wei Ding (Department of Plant Pathology and Microbiology, University of California, Riverside, USA) and Prof. Yi-Jun Qi (Tsinghua-Peking Center for Life Sciences, and School of Life Sciences, Tsinghua University, China). B. napus plants were grown in growth cabinets at 25 • C under a 16/8 h light/dark photoperiod, while Arabidopsis plants of the wildtype and mutants of RNA-silencing related genes were grown at 23 • C with a 12/12 h day/night photoperiod. Fresh sclerotia of S. sclerotiorum were cultured at 23 • C on potato dextrose agar medium (PDA) to produce mycelia, which were transferred to new PDA plates and grown for 2 days. The PDA plugs containing the advancing edge of S. sclerotiorum mycelia were removed to inoculate the plant leaves. For gene expression analyses, leaves were collected at 0, 8, and 16 h post inoculation (hpi) and frozen immediately in liquid nitrogen. Diameter of disease lesions was measured at 24 hpi and statistically analyzed using SPSS (verson19.0) by Student's t-test (p < 0.05). For disease resistance evaluation, at least 10 plants for each genotype were examined and the experiments were conducted three times independently.

Real-Time Quantitative PCR
Total RNA was extracted using Trizol reagent (Invitrogen, CA, USA) following the manufacturer's procedure. Real-time quantitative PCR (RT-qPCR) was carried out using SYBR Premix Ex Taq (TakaRa, China) on StepOne Real-Time PCR System (ABI, USA).The RT-qPCR analyses were conducted three times, with three replicates for each gene and the relative fold changes were calculated using the 2 − Ct method as described . A B. napus elongation factor gene was used as the reference gene and primers used for RT-qPCR are listed in Table S2. Significance of the differences between mean values was determined with Student's t-test (p < 0.05).

Genome-Wide Identification of DCL, AGO, and RDR Genes in B. napus
A BLASTp search was conducted against B. napus genome in NCBI database using well-characterized Arabidopsis AGO, DCL, and RDR protein sequences as query sequences. The retrieved sequences were further analyzed for domain composition. Finally, 8 DCL, 27 AGO, and 16 RDR genes were identified in B. napus genome ( Table 1). Compared with A. thaliana which contains 4 DCL, 10 AGO, and 6 RDR genes (Table S1); the members for each gene family in B. napus expanded by two times or more (Table 1). To compare the composition of these RNA silencing machinery genes between B. napus and its progenitor species B. rapa and B. oleracea, similar BLASTp and domain identification analyses were performed for these two genomes. The results showed that B. rapa genome contained 4 DCL, 13 AGO, and 6 RDR genes while B. oleracea genome carried 4 DCL, 14 AGO, and 7 RDR genes (Table S1). Comparison analysis  indicated that B. napus genome possessed all copies of DCL, AGO, and RDR genes from the two progenitor species and contained three extra copies of RDR genes ( Table 1 and Table S1).

Classification of B. napus DCL, AGO, and RDR Genes Based on Phylogenetic Analysis
In order to examine the phylogenetic relationship of DCL, AGO, and RDR families, we constructed unrooted phylogenetic trees of all BnDCL, BnAGO, and BnRDR protein sequences along with their A. thaliana, B. rapa, and B. oleracea homologs (Figure 1). The 8 BnDCLs were obviously divided into four groups as reported for AtDCLs. Each group comprised two members with one contributed by each subgenome (A and C; Figure 1A). Coincidently, B. rapa and B. oleracea contained 4 DCLs with one member for each group ( Figure 1A). According to the phylogenetic relationship and sequence homology with AtDCLs, the 8 BnDCLs were named as BnDCL1A, BnDCL1C, BnDCL2A, BnDCL2C, BnDCL3A, BnDCL3C, BnDCL4A, and BnDCL4C in accordance with their genomic localization (in A or C). Besides, BnDCLs showed high sequence similarities (from 78.7 to 87.6% for each group) to their Arabidopsis counterparts (Table S3). Based on the NJ phylogenetic tree and the protein sequence homologies with AtAGOs, the BnAGOs family consisted of 4 AGO1s (BnAGO1A1, BnAGO1A2, BnAGO1C1, and BnAGO1C2), 2 AGO2s (BnAGO2A and BnAGO2C), 2 AGO3s (BnAGO3A and BnAGO3C), 4 AGO4s (BnAGO4A1, BnAGO4A2, BnAGO4C1, and BnAGO4C2), 2 AGO5s (BnAGO5A and BnAGO5C), 2 AGO6s (BnAGO6A and BnAGO6C), 3 AGO7s (BnAGO7A, BnAGO7C1, and BnAGO7C2), 1 AGO8 (BnAGO8A), 5 AGO9s (BnAGO9A1, BnAGO9A2, BnAGO9C1, BnAGO9C2, and BnAGO9C3) and 2 AGO10s (BnAGO10A and BnAGO10C). The total B. napus AGO copies in each subgroup were corresponding to that in the two progenitors B. rapa and B. oleracea except one extra AGO9 (BnAGO9C3) and one less AGO8 ( Figure 1B). Notably, an uneven number of AGO gene copies from these three Brassica species was observed. Genomes of B. rapa and B. oleracea comprised two copies of AGO1, AGO4, and AGO9 genes, which was identical to the subgenomes A and C of B. napus except that the subgenome C of B. napus contained an extra copy of AGO9 (BnAGO9C3). In addition, B. oleracea genome and B. napus subgenome C possessed two copies of AGO7 genes. Thus, copy numbers of these AGO genes were higher in these Brassica species than in A. thaliana. Instead, B. napus subgenome C did not contained any AGO8 gene. Besides the exceptions herein described, the number (only one) of the remaining AGOs in genomes of B. rapa and B. oleracea and subgenomes A and C of B. napus was identical to A. thaliana. On the other hand, the distribution of gene members of AGO groups was generally even in the A and C subgenomes of B. napus except AGO7, AGO8, and AGO9 ( Figure 1B and Table 1). Additionally, sequence similarity between BnAGOs and Arabidopsis homologs was generally high, ranging from 60.8 to 92.7%, while the similarity among sequences of the same group of BnAGOs was 87.9-91.1% for AGO1, 86.7% for AGO2, 90.0% for AGO3, 96.8-99.6% for AGO4, 94.5% for AGO5, 98.6% for AGO6, 76.5-91.0% for AGO7, 72.4-97.3% for AGO9 and 96.1% for AGO10 (Table S3).
Collectively, genome of the allopolyploid species B. napus possesses almost all copies of DCL, AGO, and RDR genes from its progenitor species B. rapa and the RDR5 group in B. napus appears to have undergone further expansion through duplication during evolution. Furthermore, compared with Arabidopsis, some AGO and RDR genes such as AGO1, AGO4, AGO9, and RDR5 have significantly expanded in these Brassica species.

Physico-Chemical Properties and Domain Composition of B. napus DCL, AGO, and RDR Proteins
Physico-chemical properties and domain composition of B. napus DCL, AGO, and RDR proteins were generally similar to their Arabidopsis counterparts, though some differences were noticed. For DCLs, seven out of eight BnDCL proteins contained DEAD, Helicase-C, Duf283, PAZ, and RNaseIII domains as reported for Arabidopsis DCLs. The remaining one (BnDCL3C) lacked the DEAD domains, which is distinguishable from AtDCL3 in this regard (Figure 2A). Further, comparison of the DEAD domain region of all Arabidopsis and three Brassica species revealed two deletions of 31 and 9 amino acids (aa), respectively, in BnDCL3C ( Figure S1A). The correct sequence of BnDCL3C awaits further experimental confirmation. The gene length was highly similar within groups but considerably varied between groups of BnDCLs. BnDCL1s were the largest (5439 and 5430 bp) followed by BnDCL4s (4929 and 4926 bp) and BnDCL3s (4572 and 4596 bp), while BnDCL2s were the smallest (4167 and 4170 bp; Table 1). This group-wise variation in gene length is also observed in AtDCLs (Table S1).
The domain composition of BnAGOs was identical to that of AtAGOs. All BnAGOs possessed DUF1785, PAZ, and PIWI domains. Besides, all BnAGO1 proteins contained an additional Gly-rich Ago1 domain ( Figure 2B). Furthermore, we examined BnAGOs for presence of the four key active residues (DDH/H) in the PIWI domain that are responsible for the endonuclease property of AGO proteins involved in RNAi. The PIWI domain sequence alignment revealed that all 11 members of groups BnAGO1s (4), BnAGO5s (2), BnAGO7s (3), and BnAGO10s (2) possessed the four key active residues, suggesting that they might have endonuclease activity (Figure 3). The remaining BnAGOs belonging to AGO2, AGO3, AGO4, AGO6, AGO8, and AGO9 groups contained substitutions of the two H residues, which is similar to their Arabidopsis AGO counterparts. Furthermore, the length of the BnAGOs CDS varied from 2604 bp for BnAGO6 to 3261 bp for BnAGO1A2, potentially encoding 867 and 1086 amino acids, respectively ( Table 1). All BnAGOs encode for ∼100 kDa basic proteins with a pI ranging from 8.82 to 9.66. The physico-chemical properties of BnAGOs were generally similar to AtAGOs and conserved within all AGO groups.
All the 16 BnRDR proteins carried an RdRP domain ( Figure 2C). Additionally, members of BnRDR1, BnRDR2, and BnRDR6 groups bore a RBD domain as observed in the same groups of AtRDRs ( Figure 2C). Furthermore, we examined BnRDRs for presence of the catalytic motif in the RdRP domain. The sequence alignment demonstrated that all members of BnRDR1, BnRDR2, and BnRDR6 groups shared a common DLDGD motif in the catalytic domain, while all members of BnRDR3 and five out of six BnRDR5 proteins contained DFDGD motif. The exception was BnRDR5A3, which did not contain Frontiers in Plant Science | www.frontiersin.org DFD in the characteristic catalytic motif, which requires further experimental confirmation ( Figure S1B). These data indicated that BnRDR1, BnRDR2, and BnRDR6 proteins belong to RDRα class, while the BnRDR3 and BnRDR5 proteins are members of RDR γ class. This is similar to what has been observed for AtRDRs. Additionally, the length of BnRDR proteins varied from 861 to 1198 aa ( Table 1). Different groups of RDRs exhibited diverse pI-values; the pI-value of BnRDR3s and BnRDR4 was higher than 8.0, while that of the majority of the other BnRDRs was lower than 7.0 ( Table 1).
Comparison of the protein pI-value indicated that AGOs are obviously basic proteins with a pI higher than 8.8, DCLs are generally acidic proteins with a pI lower than 6.8 except BnDCL2C, while RDRs may be acidic or basic with a pI ranging from 5.7 to 8.6 ( Table 1).

Exon-Intron Organization of B. napus DCL, AGO, and RDR Genes
The exon-intron structure of DCL, AGO, and RDR genes was examined to gain more insights into their possible structural evolution. Our results for all three gene families showed that intron number was generally conserved in members of the same groups while varied significantly in different groups of the same family. For BnDCL genes, the intron number varied from 17 to 24 and BnDCL3 and BnDCL4 groups contained 3∼7 more introns than the remaining two groups (Table 1 and Figure  S2A). In the case of BnAGOs, considerable variations in intron number were observed. BnAGO groups 1, 4, 5, 6, 8, 9, and 10 comprised similar number of introns ranging from 16 to 23, while BnAGO groups 2, 3, and 7 possessed only 1∼5 introns, which were dramatically less than other BnAGO groups (Table 1 and Figure S2B). Similar to BnAGOs, remarkable intron number variation also occurred in BnRDRs. BnRDR1, BnRDR2, and BnRDR6 groups contained as few as 1∼3 introns, while BnRDR3, BnRDR4, and BnRDR5 groups carried as many as 16∼19 introns (Table 1 and Figure S2C). These group-dependent exon-intron structures were conserved for genes from both B. napus and Arabidopsis. The observation that exon-intron structure of DCL, AGO, and RDR genes was highly similar within members of the same groups but significantly divergent across different groups of the same family suggests that these gene families especially AGO and RDR families have undergone frequent gene duplication and recombination throughout evolution.

Prediction of Cis-Acting Elements in Promoter of B. napus DCL, AGO, and RDR Genes
The 1.5 kb sequences upstream of the translation initiation codon of BnDCL, BnAGO, and BnRDR genes were retrieved and searched for cis-acting elements using PLACE database. This revealed the presence of various cis-acting elements related to phytohormone response, abiotic stress and defense response in BnDCL, BnAGO, and BnRDR genes ( Table S4). All of these genes contained at least three cis-acting elements directly related to defense response, including ASF1MOTIFCAMV (S000024), GT1CONSENSUS (S000198), SEBFCONSSTPR10A (S000391), WBOXATNPR1 (S000390), and WRKY71OS (S000447). Interestingly, 21 out of these 51 genes possessed a 6-bp CGCG element (A/C/G) CGCG (C/G/T; Figure 4 and Table  S4), with BnAGO2A, BnAGO5C, BnAGO6A, and BnAGO8A containing two copies of these cis-elements. It has been revealed that calmodulin-binding transcription activators (CAMTAs) contribute to plant defense responses by binding to CGCG cis-elements of the target gene promoters and thereby regulating their expression (Yang and Poovaiah, 2002;Du et al., 2009;Nie et al., 2012). Thus, intriguingly, our result predicts that RNA silencing might be regulated by CAMTAs.

Expression Analyses Implied That CAMTA3 Might Mediate Regulation of B. napus DCL, AGO, and RDR Gene Expression in Response to Pathogen Inoculation
Owing to the high cDNA sequence identity in the same gene family, design of gene-specific primers for many of B. napus DCL, AGO, and RDR genes is unfeasible. Therefore, we chose all gene members of BnAGO4, BnRDR1, and BnDCL1 groups, for which ideal gene-specific primers could be designed, as representative to examine the expression pattern of these predicted genes in B. napus and their response to S. sclerotiorum inoculation. The analysis demonstrated that these genes were differentially expressed in leaves under normal growth conditions ( Figure 5A and Figure S3). Six of them including four BnAGO4s, BnRDR1A, and BnDCL1A, exhibited moderate to high level expression in normal leaves, while the expression of BnRDR1C1, BnRDR1C2, and BnDCL1C was only detected when a second round PCR was conducted using the products of the first round PCR as templates ( Figure 5A and Figure S3). These data demonstrated the difference in constitutive expression of these nine BnAGO4, BnRDR1, and BnDCL1 genes in B. napus leaves. Real-time quantitative PCR (RT-qPCR) was further used to gain insight into expression patterns of these genes in response to S. sclerotiorum inoculation at 0, 8, and 16 hpi. Interestingly, all of these genes, except BnRDR1A, BnRDR1C1, and BnDCL1A, were downregulated to various degrees by S. sclerotiorum inoculation. Notably, 4 AGO4s were all markedly reduced by 2.2∼10.2 folds at 16 hpi ( Figure 5B). This result suggested the possible involvement of these genes in the interactions between B. napus and S. sclerotiorum.
All these genes with the exception of 2 BnRDR1 (BnRDR1A and BnRDR1C1), carried CGCG elements in the DNA sequences upstream of their coding regions (Figure 4) and we conjectured that the lowered expression of these genes might have resulted from up-regulation of CAMTA genes in response to S. sclerotiorum inoculation. To test this hypothesis, we examined the expression of BnCAMTA3 genes after inoculation with this pathogen. BnCAMTA3 genes were selected since AtCAMTA3 has been reported to function in regulating plant defense (Du et al., 2009;Nie et al., 2012;Rahman et al., 2016). Four BnCAMTA3 genes were identified through BLASTp searches against B. napus genome database using AtCAMTA3 as query. RT-qPCR assays demonstrated that expression of two out of four BnCAMTA3 genes (BnCAMTA3a and BnCAMTA3b) increased  BnCAMTA3 genes (C) in response to Sclerotinia sclerotiorum inoculation in B. napus leaves. Expression of these genes at 8 and 16 hpi was shown relatively to that at 0 hpi. The experiments were conducted three times, each containing three replicates for all genes. The expression data were statistically analyzed using SPSS software. Significance of the differences between mean values was determined with Student's t-test. Error bars represent SD, while asterisks (*) indicate significant difference at p < 0.05. strongly by 14.7 and 10.4 folds, respectively, at 8 h after S. sclerotiorum inoculation (Figure 5C). Given that CAMTA3 negatively regulates plant disease resistance through direct binding and subsequent repression of target defense-related genes (Du et al., 2009;Nie et al., 2012), our results indicate that RNA silencing might be regulated by CAMTAs during B. napus-S. sclerotiorum interactions.
Arabidopsis dcl, ago, and rdr Mutants Commonly Exhibited Altered Susceptibility to S. sclerotiorum The observation in this study that expression of BnDCL1, BnAGO4, and BnRDR1 genes altered significantly in response to S. sclerotiorum inoculation and our previous findings that Arabidopsis ago1 and ago2 mutant plants exhibit enhanced susceptibility to S. sclerotiorum (Cao et al., 2016) prompted us to assess the possible role of DCL, AGO, and RDR-mediated RNA silencing in plant resistance to this pathogen. A total of 13 A. thaliana dcl, ago, and rdr mutants were tested. They included three dcl, seven ago, and three rdr mutants, all in Col-0 background except dcl1-9 and ago4-1 which are derived from Ler ecotype. Based on the phenotypes after S. sclerotiorum inoculation, the mutants could be divided into three classes. Class I comprised those exhibiting enhanced susceptibility, Class II contained mutants displaying enhanced resistance, while Class III included mutants showing similar disease phenotypes to wildtype plants (Figure 6). The mutants dcl4-2, ago9-1, rdr1-1, rdr6-11, and rdr6-15 were more susceptible to S. sclerotiorum challenge by showing more severe necrosis and larger size of disease lesions FIGURE 6 | Susceptibility to S. sclerotiorum in Arabidopsis dcl, ago, and rdr mutants. (A) Disease symptoms of dcl, ago, and rdr mutants and wild type plants after inoculation with S. sclerotiorum. Photographs were taken at 24 hpi. I, mutants exhibiting enhanced susceptibility; II, mutants displaying enhanced resistance; and III, mutants showing similar disease phenotypes to wild-type plants. (B) Statistical analysis of disease severity. The inoculation analysis was performed three times, each in at least 10 plants for all mutants. The lesion size data were statistically analyzed using SPSS software. Significance of the differences between mean values was determined with Student's t-test. Error bars indicate SD, while asterisks (*) indicate significant difference at p < 0.05. than wild-type plants (Figure 6). By contrast, the mutant dcl1-9 was found to be more resistant since it displayed less necrosis and smaller size of disease lesions when compared with Ler plants upon S. sclerotiorum inoculation (Figure 6). However, the dcl2-1, ago3-1, ago4-1, ago5-1, ago6-1, ago7-1, and ago10-3 mutant plants exhibited similar disease susceptibility as the wild-type plants with respect to severity of necrosis and size of disease lesions (Figure 6). These data demonstrate the involvement of DCL, AGO, and RDR-mediated RNA silencing in plant resistance to the necrotrophic fungal pathogen S. sclerotiorum.

DISCUSSION
RNA silencing is a versatile molecular mechanism to regulate diverse biological processes through altering transcript accumulation of genes essential to these processes. In this study, we identified three gene families encoding key components of RNA silencing machinery in B. napus and its progenitors B. rapa and B. oleracea. We demonstrated the important role of RNA silencing in plant resistance to the fungal pathogen S. sclerotiorum and indicated a possible new molecular mechanism to regulate RNA silencing in response to this pathogen. Our results provide insights into composition, function, and mechanism of RNA silencing in plants.

RNA Silencing Machinery in B. napus
DCL, AGO, and RDR are key components of plant RNA silencing machinery. Our results demonstrate that B. napus possesses 8 BnDCL, 27 BnAGO, and 16 BnRDR genes. This gene copy numbers are much larger than those in A. thaliana although the two species belong to the same family (Brassicaceae). However, considering that B. napus is tetraploid while A. thaliana is diploid, it is clear that each haploidy in the two species carries the same number (4) of DCL genes but different number of AGO and RDR genes. The haploid B. napus contains more AGO and RDR genes than A. thaliana. B. napus is an allotetraploid from crossing between B. oleracea and B. rapa, followed by chromosome doubling (Chalhoub et al., 2014). Our results demonstrate that B. napus genome includes almost all copies of these RNA silencing machinery genes from its progenitor species B. rapa and B. oleracea. In addition, B. napus genome contains three extra copies of RDR5 genes, indicating that the RDR5 group in B. napus appears to have undergone further expansion through duplication during evolution. Moreover, compared with A. thaliana, some AGO and RDR genes such as AGO1, AGO4, AGO9, and RDR5 have significantly expanded in these Brassica species. The B. napus subgenomes A and C bear more RDR5 genes but similar AGO genes compared with genomes of its progenitor species B. rapa and B. oleracea, suggesting that the expansion of RDR5 genes is likely to occur after B. napus formation while that of AGO1, AGO4, and AGO9 may have occurred before B. napus formation. Additionally, members of AGO7, AGO8, and AGO9 as well as RDR4 are unevenly distributed in the A and C subgenomes of B. napus. This asymmetric distribution may have resulted from homeologous exchanges, which is also most likely the source of further expansion of these two families and it is consistent with the previous report (Chalhoub et al., 2014). Considering that B. napus genome contains at least two copies of each DCL, AGO, and RDR genes except AGO8, and that all four members of BnAGO4 are expressed constitutively and/or in response to pathogen inoculation ( Figure 5) and are thus most probably functional, it is interesting to understand whether these genes functions redundantly or differentially.

Function of RNA Silencing in Plant Resistance to Fungal Pathogens
It is well-known that RNA silencing can protect plants from viral infection (Ding and Voinnet, 2007;Llave, 2010;Incarbone and Dunoyer, 2013;Wu et al., 2015). This antiviral immunity involves production of virus-derived small interfering RNAs (viRNAs) and results in specific silencing of viruses by viRNAguided effector complexes. Apart from defense against viruses, RNA silencing also plays a role in plant defense against bacterial pathogens (Katiyar-Agarwal et al., 2006;Navarro et al., 2006;Voinnet, 2008;Robert-Seilaniantz et al., 2011). Similar to viruses, bacteria have also developed mechanisms to suppress RNA silencing in order to infect successfully (Navarro et al., 2008). More recently, functions of RNA silencing in plant resistance to fungal pathogens are being revealed. Plant miRNAs are differentially expressed in response to inoculation with fungal pathogens, such as Erysiphe graminis (Xin et al., 2010), Fusarium virguliforme (Radwan et al., 2011), V. dahliae (Yin et al., 2012;Yang et al., 2013), V. longisporum (Shen et al., 2014), M. oryzae (Li et al., 2014), and B. cinerea (Jin and Wu, 2015). More importantly, mutants of key components of the RNA silencing machinery exhibit altered susceptibility to fungal pathogens including two species of Verticillium (Ellendorff et al., 2009;Shen et al., 2014). Moreover, it has been reported that B. cinerea produces small RNAs to suppress plant defense by hijacking host RNA interference pathways (Weiberg et al., 2013). It is interesting to elucidate the function and mechanism of RNA silencing in the interactions between the important oil crop B. napus and the devastating fungal pathogen S. sclerotiorum. Previously, we have identified the miRNAs involved in this plant-pathogen interaction, many of which targeting genes involved in plant defense. Besides, three miRNAs (ath-miR168a_1ss21AC, aly-miR403a-3p_L+1, and bna-miR403) target AGO1 and AGO2 which are two key components of RNA silencing. We further found that Arabidopsis ago1-27, ago1-33, and ago2-1 mutant plants exhibit enhanced susceptibility to S. sclerotiorum (Cao et al., 2016), thus providing a clue to the important roles of RNA silencing in the interactions between B. napus and S. sclerotiorum. In this study, we identified gene families encoding DCL, AGO, and RDR, three key components of RNA silencing in B. napus. We found that these genes, represented by all members of three groups of genes (4 BnAGO4, 3 BnRDR1, and 2 BnDCL1 genes), are differentially expressed in response to S. sclerotiorum inoculation ( Figure 5B). Our results showed that the expression divergence was present among members belonging to the same gene group after S. sclerotiorum inoculation (Figure 5B), which suggested that gene members within the same groups may have diverse functions in this plant-pathogen interaction. AtAGO4 is one of the critical components in the transcriptional genesilencing pathway associated with siRNA that directs DNA methylation Qi et al., 2006) and is required for resistance to P. syringae (Agorio and Vera, 2007). AtRDR1 is elicited by salicylic acid (SA) treatment and viral infection (Yu et al., 2003;Diaz-Pendon et al., 2007;Qi et al., 2009). Together, these results suggest that AGO4 and RDR1 may also contribute to B. napus defense against S. sclerotiorum. Furthermore, we assessed susceptibility toward S. sclerotiorum in 13 dcl, ago, and rdr A. thaliana mutants. As many as six dcl, ago, and rdr mutants exhibit altered susceptibility (Figure 6). These mutants include dcl4-2, ago9-1, rdr1-1, rdr6-11, rdr6-15, and dcl1-9. The mutated genes encode different RNA-silencing components. This is similar to what was found in the study involving another fungal pathogen V. dahliae (Ellendorff et al., 2009). Collectively, these results indicate that RNA silencing may be involved in interactions between plants and fungal pathogens. However, how these DCLs, AGOs, and RDRs regulate resistance to S. sclerotiorum is unclear. Recently, it has been reported that rice AGO18 functions through regulating AGO1 accumulation by direct binding with the AGO1-targeted miRNA (miR168; Wu et al., 2015). Whether a similar mechanism exists for AGOs or even RDRs and DCLs in Arabidopsis and B. napus is worthy of further study.

Regulation of RNA Silencing Machinery
A remaining challenge for RNA silencing is to dissect the mechanism through which the RNA silencing pathway itself can be regulated. Up to now, little of this mechanism has been uncovered except the evidence that the expression of components of the RNA silencing pathway is subject to negative feedback regulation by their own miRNA products. For example, miR162 targets DCL1, miR168 targets AGO1, and miR403 targets AGO2. However, there is no more information regarding regulatory mechanisms for the other RNA silencing components. Considering the importance of RNA silencing pathway, we are curious to know whether other mechanisms exist that confer additional layers of regulation on the RNA silencing machinery. Intriguingly, in present study, we found that substantial number of B. napus DCL, AGO, and RDR genes (21 out of a total of 51) possessed CAMTA/SR binding sites [(A/C/G) CGCG (C/G/T)] in their promoter sequences (Figure 4 and Table  S4). CAMTAs, especially CAMTA3, contribute to plant defense responses by direct binding and thereby regulating the expression of the target genes (Yang and Poovaiah, 2002;Du et al., 2009;Nie et al., 2012). Therefore, we suspect that the expression of these CGCG-element-containing RNA silencing components may be regulated by CAMTAs. To confirm this possibility, we examined the expression of CAMTA genes and these CGCGelement-containing RNA silencing genes in B. napus in response to S. sclerotiorum inoculation. The results of this expression analysis indicate that S. sclerotiorum inoculation strongly induced the expression of BnCAMTA3 genes while it significantly suppressed that of many CGCG-element-containing BnAGO, BnDCL and BnRDR genes ( Figure 5). Moreover, another work in our laboratory has revealed that Atcamta3 mutant plants exhibit enhanced resistance to S. sclerotiorum (Rahman et al., unpublished data). Taken together, our results suggest that RNA silencing might be regulated by CAMTA3. Nevertheless, further confirmation of CAMTA binding activity with CGCGelement-containing RNA silencing genes by other assays such as EMSA and its effect on expression of these genes will provide more straightforward evidence to support this intriguing likely mechanism for regulation of RNA silencing machinery.
It is well-known that the functions of CAMTAs are dependent on their interaction with Ca 2+ /CaM (Choi et al., 2005;Du et al., 2009) and these genes, especially CAMTA3, are widely involved in plant defense (Yang and Poovaiah, 2002;Du et al., 2009;Nie et al., 2012;Rahman et al., 2016). For instance, knockout of AtCAMTA3 leads to increased accumulation of salicylic acid and enhanced host disease resistance to both bacterial (Du et al., 2009) and fungal pathogens (Nie et al., 2012) and nonhost resistance to bacterial pathogen (Rahman et al., 2016) but reduced resistance against insect herbivores (Qiu et al., 2012). Similarly, one rice CAMTA mutant (oscbt-1) exhibits enhanced resistance to blast fungal pathogen and leaf blight bacterial pathogen (Koo et al., 2009). Given that expression of some RNA silencing machinery may be regulated by CAMATs and both of CAMTA as well as RNA silencing contribute to plant resistance, it will be intriguing to explore whether the role of CAMTAs in regulating plant defense response can be attributed to its regulation of RNA silencing and thus opening a possibility to link Ca 2+ signaling and RNA silencing together to provide a novel facet of Ca 2+ signaling in regulation of plant disease resistance.

CONCLUSIONS
The B. napus genome possessed 8 Dicer-like (DCL), 27 Argonaute and 16 RNA-dependent RNA polymerase (RDR) genes. The B. napus genome expanded RDR5 genes compared with that of its progenitors B. oleracea and B. rapa and all genomes of three Brassica species expanded AGO1, AGO4, and AGO9 genes compared with the Arabidopsis genome. B. napus DCL, AGO, and RDR genes widely (21 out of 51) harbored a CAMTA-binding site (CGCG box) in their promoter regions. Inoculation with S. sclerotiorum strongly induced the expression of BnCAMTA3 genes while significantly reduced that of many CGCG-containing RNA silencing component genes. Our results suggested that RNA silencing machinery might be targeted by CAMTA3. Furthermore, mutant analyses demonstrated that six out of 13 dcl, ago, and rdr mutants exhibited altered susceptibility to S. sclerotiorum. These results indicate the important role of RNA silencing in plant resistance to this devastating fungal pathogen.

AUTHOR NOTE
During the review process of this manuscript, a paper on bioinformatics identification of DCL, AGO, and RDR families in Brassica napus was accepted for publication .

AUTHOR CONTRIBUTIONS
The project was coordinated by X-ZC. J-YC, and Y-PX conducted the bioinformatics and phylogenetic analyses. J-YC, Y-PX, and HR carried out the gene expression assays. J-YC, WL, and S-SL performed disease resistance evaluation analyses. J-YC designed and performed the statistical analysis. X-ZC conceived of the study, and participated in its design and coordination. X-ZC and J-YC prepared the manuscript. All authors read and approved the final manuscript.