Identification and Analysis of miRNAs and siRNAs in Botrytis cinerea

Small RNAs in Botrytis cinerea were analyzed via high-throughput sequencing on BGISEQ-500 platform. A total of 8 novel miRNAs and 110 novel siRNAs were predicted. Sequence information, construction, length distribution, base bias and expression levels of miRNAs and siRNAs were determined as well. Through GO and KEGG enrichment analysis, the miRNA target genes are mostly located in membrane and organelle, possessed binding and catalytic activities, and involved in signal transduction and carbohydrate metabolism. The results will provide a theoretical foundation for understanding the developmental and pathogenic mechanisms of B. cinerea at the transcriptional level.

morphological observation and rDNA-internal transcribed spacer analysis in the preliminary experiments. The fungus was cultivated on a potato dextrose agar (PDA) plate at 25˚C. To evaluate the mycelial growth, a mycelial agar disk (0.5 cm in diameter) was placed in the center of a 9 cm diameter petri dish containing 25 mL PDA. During 10 days of culturing at 25˚C, radial growth was measured daily by the decussation method. Meanwhile, the mycelial phenotype was microscopically observed using a Nikon Eclipse Ni-U microscope (Nikon, Japan). Tomato fruit at commercial maturity were bought from local market. After disinfection by 2% sodium hypochlorite, fruit were wounded at the equator utilizing a sterile needle and inoculated with 10 μL spore suspension (1 × 10 4 spores/mL). The lesion diameter of inoculated fruit was recorded in every two days.

Small RNA Sequencing
Fresh spore suspensions were prepared by flooding the sporulating cultures of B. cinerea with sterile distilled water containing 0.05% (v/v) Tween-20. A suitable aliquot of the suspension was added to 100 mL potato dextrose broth (PDB) with a final concentration of 1 × 10 6 spores/mL. After 24 h of shaking culture at 25˚C, the harvest spores and mycelia were washed twice with sterile distilled water, and quickly frozen in liquid nitrogen. A biological repeat was carried out, and the two samples were mixed with the same weight. Small RNA sequencing was performed using BGISEQ-500 technology of Beijing Genomics Institute (BGI) Co., Ltd. Briefly, small RNA was separated from total RNA by PAGE gel, and linked with a 5'-adenylated, 3'-blocked single-stranded DNA adapter at 3' end. After RT primer hybridization, the 5' adaptor was linked and the first strand cDNA was synthesized. To enrich cDNA of 100 -120 bp, PCR amplification and PAGE gel separation were carried out. Finally, the library was quantified and pooling cyclization was performed. Each step was under strict quality control.

Bioinformatic Analysis
After sequencing, the original data were filtered by removing low-quality reads, adaptors and other contaminants to get clean reads. The remaining tags (clean tags or clean reads) were stored in FASTQ format. Bowtie2 was used to map clean reads to the reference genome and to other sRNA [12], and Cmsearch was used for Rfam [13]. sRNA classification followed the priority rule as MiRbase > pirnabank > snoRNA > Rfam > other sRNA. Novel miRNA was predicted by exploring the characteristic hairpin structure of miRNA precursor in miRDeep2 and miRA [14] [15]. siRNA predicting followed the criteria that siRNA was a 22 -24 nt double-stranded RNA, each strand of which was 2 nt longer than the other [16]. The small RNA expression level was calculated by using TPM (transcripts per million). TPM = C × 10 4 /N. C means miRNA counts number in a sample, and N means total reads number that mapped to the genome [17]. TAPIR  TargetFinder were used to predict miRNA target genes [18]. The default parameters are as follows, TAPIR: --score 5 --mfe_ratio 0.6; TargetFinder: -c 4.
WEGO software was used to perform Gene Ontology (GO) enrichment analysis of miRNA target genes. KEGG database was used to perform pathway enrichment analysis of miRNA target genes. A scatter plot of the KEGG enrichment results and a bar plot of the KEGG terms were generated. A corrected P value ≤ 0.05 was taken as a threshold.

Statistical Analysis
Data were pooled across independent repeat experiments. Correlation analysis was carried out to test the significance of the relationship between two or more variables utilizing Statistical Product and Serviced Solutions software (SPSS, USA). Analysis of variance (ANOVA) was used to compare more than two means. Mean separations were analyzed using Duncan's multiple range test.
Differences at P < 0.05 were considered to be significant.

Phenotype and Activity of B. cinerea
In a previous experiment, the isolate has been identified as Botrytis cinerea based on morphological and genetic characteristcs. Typically, sclerotia of B. cinerea commence growth to produce conidiophores and multinucleate conidia, serving as a primary source of inoculum. Microconidia are also formed as spermatia. The sexual cycle involves the spermatization of sclerotia, leading to the production of apothecia and asci with eight binucleate ascospores. Conidia generation follows a cycle of initiation, production and dissemination that is controlled by changes in light, humidity and temperature. Dry-inoculated conidia produced one or two short germ tubes and no obvious terminal appressoria, whereas conidial within suspension possessed much longer germ tube and extensive secretion of an extracellular matrix. On PDA plate, B. cinerea initially formed white colony turning gray or grayish brown ( Figure 1(A)). The mycelial extension was fast and the colony could cover the petri dish in 8 days ( Figure   1(C)). The hyphae were septate and intertwined (Figure 1(B)). Conidiophores were gray brown, ellipsoidal or ovate.
B. cinerea is responsible for a very wide range of symptoms. Soft rots, accompanied by collapse and water-soaking of parenchyma tissues, followed by a rapid appearance of grey masses of conidia are the most typical symptoms on fruit.
Soft rotting of mature tomato fruit occurs mainly in postharvest stage. In the present study, at the temperature of 25˚C with high relative humidity, the infection of B. cinerea ontomato fruit was rapid, resulting in softening of the flesh and a browning skin. On the fruit surface, massive mycelia were visible. Within about a week, the fruit was completely rotten (Figure 1(D)). Therefore, the fungus used in this study presented a typical phenotype of B. cinerea, and exhibited a normal growth and pathogenicity.

Small RNA Sequencing Results
The small RNAs were sequenced using BGISEQ-500 technology. A total of 27,872,641 raw tags containing 1,006,245 low quality tags, 875,820 invalid adapter tags, 774 polyA tags, 1,934,729 short valid length tags and 24,055,073 clean tags. The base percentage composition of clean tags was shown in Figure 2(A).
The percentage of clean tags, of which quality was more than 20, was 98.80% of all clean tags (Figure 2 The proportion of all kinds of sRNA was shown in Figure 3(A) and the genome distribution of tags was shown in Figure 3(B).

Prediction of miRNAs and siRNAs
Only miR-466i-5p as known miRNA was found in B. cinerea. It belonged to mmu-miR-466i-5p family and its hairpin contained 300 bases. Eight novel miRNAs named as Bc-mir1 to Bc-mir8 were predicted. Their expression levels  were significantly different. The expression of Bc-mir1 was about 100 times higher than that of Bc-mir8. The sequence information of novel miRNAs was shown in Table 1. The stem loop structure of precursors was shown in   indicate the predicted miRNAs Bc-mir1 to Bc-mir8.
composition of predicted miRNAs and siRNAs was shown in Figure 5. For novel miRNAs, the first base distribution was G, C, A and U in number from high to low.
There were no obvious rules for the base distribution due to the small number. The length of novel siRNAs was mostly less than 21 nt. The first base distribution was U,  A. In addition, the twenty-third base of 14 siRNAs did not contain base A.

Enrichment Analysis of miRNA Target Genes
Software TAPIR and TargetFinder were used to find the target genes of miRNAs in B. cinerea. The intersection results of two softwares contained 300 target genes, and there were 642 target genes in union results (Table S2). The GO enrichment result was shown in Figure 6. Most of the target genes located in membrane and organelle, and possessed binding activity, catalytic activity or transporter activity. The target genes involved in many biological processes including cellular process, biological regulation, metabolic process, single-organism process, etc. The KEGG enrichment result was shown in Figure 7. The items with a larger number of target genes were carbohydrate metabolism, signal transduction, and transport and catabolism. In addition, statistics of pathway enrichment indicated that many target genes were related to the MAPK signaling pathway, amino sugar and nucleotide sugar metabolism, and endocytosis.
Meanwhile, target genes involved in the AGE-RAGE signaling pathway, alpha-linolenic acid metabolism, and phosphatidy linositol signaling system possessed higher rich factor value ( Figure 8).
Y. Y. Liu et al. genes. The Y axis represents GO term. All GO terms are grouped into three ontologies: blue is for biological process, green is for cellular component and red is for molecular function.

Discussion
RNAi is a highly conserved biological process in diverse species. It has in com-    and characterized two novel proteins, SAD-4 and SAD-5. Both of them were required for meiotic silencing by unpaired DNA and siRNA generation in Neurospora crassa [31]. Deshmukh and Purohit (2014) used synthetic siRNAs to down-regulate secondary metabolite genes hmgR and fpps, which led to bikaverin overproduction in Fusarium sp. HKF15 [32]. Yu et al. (2014) showed that convergent antisense transcription and availability of the Dicer ribonuclease were key determinants for heterochromatin formation and catalyzing siRNA generation in Schizosaccharomyces pombe [33]. Yu et al. (2021) found that trimethylguanosine synthase 1 (Tgs1) played critical roles in Swi6/HP1-independent siRNA production and establishment of heterochromatin in fission yeast [34].

Conclusion
In