Classification of small RNAs response to H. parasitica in non-heading Chinese cabbage
To identify miRNAs responsive to downy mildew in non-heading Chinese cabbage, 15 sRNA libraries were constructed by using leaf samples collected from five different time points (0, 6, 12, 24 or 48 hpi). High-throughput sequencing generated of mainly 18 to 30 nucleotides in size raw reads for each library of ranged between 12.5 million and 31.4 million (Table 1), which presented a deep resource for extensive discovery of small regulatory RNAs. In total, >280 million raw reads were obtained (Table 1) and processed for miRNA identification. After removing poly(A) tails, short and low-quality reads, and adaptor contamination, a total of clean reads for each library of ranged between 12.1 million and 31.1 million (Table 1). Thereafter, we mapped clean reads to known RNAs in the Silva database, GtRNAdb database, Rfam database and Repbase database (Table 2), filtering rRNA, tRNA, scRNA, snRNA, snoRNA and repetitive sequences. rRNAs were the most abundantly annotated RNAs, accounting for about half of total reads except 3 libraries of 48 hpi, and the unannotated RNAs were potentially novel regulatory small RNAs (Table 2). Then, miRBase database (v22) was utilized to identify known miRNAs. The number of unannotated clean reads that mapped to miRBase were ranged between 458535 and 4961663 (Table 2), and the number of known miRNAs were ranged between 157 and 191 (Table 1). The novel miRNAs (Table 1) were obtained after miRDeep-P [29] analysis by using the unannotated reads that unmapped to miRBase, and the number of novel miRNAs were ranged between 91 and 112 (Table 1). The formation of mature bodies is realized by the shear of Dicer/DCL enzyme and mature plant miRNAs are always 21 nt or 24 nt [30]. The 21- and 24-nt sequences were dominant in all libraries, and the 21-nt ones were the most abundant (Fig. 1a, b).
Characteristics of miRNAs response to H. parasitica in non-heading Chinese cabbage
In total, we identified 197 known miRNAs (Table S2),117 of them not described in B. rapa (Table S2), and 112 novel miRNAs (Table S3). The number of miRNAs in length of 21/24nt were 140 and 81 respectively for the known miRNAs and novel miRNAs (Fig. 2a, b; Table S2, S3), respectively. Based on the sensitivity and specificity of Dicer/DCL at recognizing and cleaving pre-miRNAs [30], miRNA bias analysis was conducted. The results indicated that the mature miRNA sequences had A-U contents greater than the G-C contents throughout all transcripts except at position 7, 15, 17, 19, 21, 22, 23 and 24 in known miRNAs, and 3, 4, 6, 8, 9, 14, 18, 19 and 21 in novel miRNAs (Fig. 3a,b). The frequencies of U appeared initially were ranged from 12.9% to 100% in known miRNAs (Fig. 4a), while they were 2.0% to 100% in novel miRNAs (Fig. 4b). Approximately 40.7% were U at the first nucleotide for 21-/24-nt miRNAs, and all first nucleotides were U for 19-nt miRNAs in known miRNAs (Fig. 4a). Approximately 21% were U at the first nucleotide for 21-/24-nt miRNAs, and all the first nucleotides were U for 18-nt miRNAs in novel miRNAs (Fig. 4b).
Differential expression of miRNAs in response to H. parasitica
To identify the response pattern of miRNAs to H. parasitica in non-heading Chinese cabbage, we compared the abundance of miRNAs in the inoculated libraries (6 hpi, 12 hpi, 24 hpi and 48 hpi) compared with the libraries at 0 hpi. Differential expressed miRNAs were identified by using DESeq2 [31] with |log2(FC)| ≥ 1 and FDR ≤ 0.05 as the significance cutoffs. The nomenclature was following the previous reports [11, 32]. A known miRNA was named as bra-miR#-5p or -3p, and the opposite strand was named as bra-miR#*-3p or -5p. A novel miRNA was named as bra-miRn#-5p or -3p, and the opposite strand was named as bra-miRn#*-3p or -5p. Total of 8 known miRNAs and 8 novel miRNAs were identified to be differentially expressed in non-heading Chinese cabbage response to H. parasitica (Table 3; Fig. 5; Table S4). To illustrate the differences in miRNAs expression levels after H. parasitica-inoculation, a heat map was constructed using Heatmap Illustrator (HemI) on the platform BMKCloud (www.biocloud.net) based on the TPM method [33]. At 6 hpi, only bra-miRn47-3p showed up-regulated expression; at 12 hpi, 6 miRNAs were down-regulated while bra-miRn38-3p was up-regulated; at 24 hpi, 4 miRNAs showed up-regulated expression and 4 miRNAs were down-regulated; at 48 hpi, 6 miRNAs were up-regulated and 2 miRNAs showed down-regulated expression (Table 3; Table S4). The differentially expressed miRNAs at different time points suggesting that these miRNAs might be induced or repressed in the response to H. parasitica. Eight of the differentially expressed miRNAs were assigned to known miRNA families (Table 3).
Prediction and function analysis of miRNA target genes
The software TargetFinder (1.6) were employed to predict the miRNA targets based on gene sequence information of known miRNAs and novel miRNAs in Brassica genomic Scaffold sequence (http://brassicadb.org/). Total of 2,722 target gene were predicted as target genes of 183 known miRNAs, and 1,662 genes were predicted as that of 85 novel miRNAs (Table 4). The entire set of miRNA target genes was annotated against the NR (ftp://ftp.ncbi.nih.gov/blast/db/), Swiss-Prot (http://www.uniprot.org/), GO (http://www.geneontology.org/), COG (http://www.ncbi.nlm.nih.gov/COG/), KEGG (http://www.genome.jp/kegg/), KOG (http://www.ncbi.nlm.nih.gov/KOG/), and Pfam databases (http://pfam.xfam.org/) using BLAST with a cutoff of E-value < 1E-05. Among the 4,384 target genes, 4,360 were annotated (Table 4, 5). 3201 target genes were annotated in GO database, of which the number of target genes at length between 300 and 1000 bp was 778, and that was 2389 at the length over 1000 bp. Only 1632 target genes annotated in KEGG database, of which 385 were at the length between 300 and 1000 bp, and 1235 at the length over 1000 bp (Table 5).
To have a better understanding of differentially expressed miRNAs, the target genes of these miRNAs were subjected to GO enrichment analysis by Blast2GO program with default parameters. The results showed that these predicted targets could be classified into 8 molecular functions, 9 biological processes and 3 cellular components (Fig. 6). Most of them, which were classified as the binding category (9 targets in oxygen binding and 11 targets in heme binding), encode transcription factors. KEGG pathway enrichment analysis was carried out as well. 30 pathways were enriched with target genes of significantly differentially expressed miRNAs, which were divided into five biological processes, including cellular processes (12.5%), environmental information processing (5.36%), genetic information processing (35.73%), metabolism (44.69%) and organismal system (1.79%) (Fig. 7). Among these pathways, miRNA surveillance pathway involved the most target genes (Bra001821, Bra028900, Bra029681, Bra040031 and Bra005265), while plant-pathogen interaction pathway only contained one target gene (Bra041056) (Fig. 7; Table S4).
qRT-PCR analysis of predicted target genes
To explore the relationship between differentially expressed miRNAs and their target genes, expression patterns of several predicted target genes were analyzed via qRT-PCR. Bra036138 were predicted as target genes of bra-miR400b-3p (Table S4). When the abundance of bra-miR400b-3p were decreased at 12 hpi and 48 hpi (Table 3), the expression level of Bra036138 were increased (Fig. 8). Bra005264 were predicted as target genes of both bra-miRn91a-3p and bra-miRn91b-3p (Table S4). When the abundance of bra-miRn91a-3p and bra-miRn91b-3p were increased at 48 hpi (Table 3), the expression level of Bra005264 were increased (Fig. 8). These suggested that miRNA-mediated silencing of their potential target genes occur during the process of defense response in non-heading Chinese cabbage after H. parasitica infection. Bra000540 and Bra024465 were predicted as target genes of bra-miRn28-5p, but they had different expression pattern (Fig. 8). When the abundance of bra-miRn28-5p was increased at 24 hpi (Table 3), the expression level of Bra024465 showed increased while Bra000540 didn’t (Fig. 8). It indicated that these target genes may not be regulated only by miRNAs upon downy mildew infection.