MicroRNAs in the response of non-heading Chinese cabbage to Hyaloperonospora parasitica (downy mildew) infection

Background MicroRNAs (miRNAs) are a subset of endogenous small non-coding ~21 nucleotide RNAs, which are playing important roles in development, abiotic stress, and pathogen responses. Downy mildew, caused by Hyaloperonospora parasitica , is a severe fungal disease in non-heading Chinese cabbage ( Brassica rapa ssp. chinensis ). However, miRNAs and their targets involved in the response of non-heading Chinese cabbage to downy mildew were still unclear. Results Here, fifteen small RNA libraries were generated from non-heading Chinese cabbage leaves at different time points after H. parasitica inoculation. By deep sequencing, we identified 197 known miRNAs and 112 novel miRNAs. According to data analysis, 8 known miRNAs and 8 novel miRNAs were found to be responsive to downy mildew infection in non-heading Chinese cabbage. Expression levels of 12 predicted miRNAs target genes were determined by reverse transcription-quantitative polymerase chain reaction (qRT-PCR). In this study, we proposed a regulatory network showing that non-heading Chinese cabbage defense response to downy mildew by repressing cell wall pectin, lipid and cellulose catabolic process ( PME12 , PME17 , PLAIVA and GDPDL4 ) and callose deposition ( CYP79B2 ), repressing indole glucosinolate metabolic process ( CYP81F3 ), repressing respiratory burst ( PUB24 ), repressing programmed cell death ( CRK5 ), and inducing defense related gene expression ( ROS1 , LECRK-VII.2 , BAK7 and WRKY48 ). Conclusions This research identified miRNAs associated with downy mildew infection in non-heading Chinese cabbage. Target genes and network analysis contributed to understanding the interaction between plant and pathogen.


Background
Brassica rapa is a species that contains multiple oil and vegetable crops [1]. Non-heading Chinese cabbage (B. rapa ssp. chinensis) is one of them and widely cultivated in Southern China [2]. Downy mildew, caused by Hyaloperonospora parasitica, is a common disease in brassicas [3][4][5][6]. Applying fungicide is an efficient way to control downy mildew [7], but it might increase the probability of pathogen population forming fungicide resistance that makes the disease management even worse [8]. Therefore, it is necessary to develop more disease-resistant varieties, which needs a better understanding of how plants response to downy mildew infection.
In plants, small RNA-guided post-transcriptional regulatory mechanisms play important roles in many aspects, including flowering regulation [9], hormone responses [10,11] and epigenetic control of transposable elements [12]. MicroRNAs (miRNAs) are small regulatory RNAs and found in diverse eukaryotic lineages [13]. They are produced from pri-miRNAs, containing a fold-back structure, via DICER-LIKE proteins processing [14]. Mature miRNAs are incorporated into an RNA-induced silencing complex (RISC) to interact with the complementary sites of the target gene transcripts, which negatively regulate the target genes expression by degrading or repressing the RNA transcripts [13,15]. Besides regulating a range of essential cellular and biological processes, miRNAs have been shown to play crucial roles in responses to abiotic and biotic stresses, such as nutritional deficiency [16], drought [17], salinity [18], cold [19], heat [20], oxidative stress [21], heavy metal stress [22], and resistance to disease [23,24]. In Arabidopsis, overexpression of miR393 enhances plant resistance to bacterium Pseudomonas syringae pv. tomato DC3000 by decreasing transcripts of TIR1, AFB2 and AFB3 that are involved in auxin signaling pathway [25]. miR398 is downregulated by Mildew resistance locus a (Mla) in barley, and thus the accumulation of SOD1 targeted by miR398 induces hypersensitive reaction after powdery mildew infection [26].
Some achievements have been made in the study of miRNAs in Brassica rapa. Yu et al. [20] showed that miR398 was downregulated and its target gene CSD1 was upregulated under heat stress. Huang et al. [11] found that miR159 participated in vernalization. miR156 and miR166 were involved in heading process in Brassica rapa [27,28]. Till now, there have been no reports on miRNAs associated with downy mildew infection in non-heading Chinese cabbage. Here, using Solexa Analyzer, we identified a set of small RNAs from downy mildew infected leaves of a resistant non-heading Chinese cabbage cultivar. Based on the small RNA sequencing and target genes prediction, we proposed a miRNA-mediated regulation model in non-heading Chinese cabbage upon H. parasitica infections.

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 upregulated; 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).
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.

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 miRNAmediated 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.

Discussion
In this study, over twelve million clean reads were obtained from each library ( Table 1) and most of them were about 20-to 24-nt long (Fig. 1a) that are common lengths of miRNAs in plant [32]. In total, 197 known miRNAs 117 of them not described in B. rapa, and 112 novel miRNAs were identified in downy mildew-infected non-heading Chinese cabbage by deep sequencing (Fig. 2; Table S2, S3).
Compared to Arabidopsis, rice or other model plants, information about miRNA in Brassica rapa is limited in miRBase (http://www.mirbase.org). This might be the reason for why so many known miRNAs not described in B. rapa were identified in our research. The number of known miRNAs were ranged between 157 and 191 ( Table 1). The novel miRNAs were obtained after miRDeep-P analysis by using the unannotated reads that unmapped to miRBase, and the number of novel miRNAs were ranged between 91 and 112 (Table 1). Compared with 0 hpi, 6 hpi, 12 hpi and 24 hpi, the number of known and novel miRNAs detected at 48 hpi was less than several other time points (Table 1). In addition, the number of miRNAs detected at 6 hpi was similar with the number of miRNAs detected at 12 hpi, but the number of differentially expressed miRNAs in response to downy mildew infection was only 1 at 6 hpi compared to 0 hpi while 7 at 12 hpi (Table 3). Since the germination time of downy mildew spores vary [34,35], we assume that the state of spores at 6 hpi and 12 hpi were different.
Two targets of bra-miR400b-3p were BAK7 and WRKY48 (Table S4), which responded to fungi invasion [43]. One target gene of bra-miRn38-3p was PLAIVA (Table S4), which was involved in the plant defense through free fatty acid and lysophospholipid signaling molecules [44]. miRNA bra-miR9552b-3p were predicted to target GDPDL4 (Table S4), which coding a protein SHV3 play important role in cell wall organization [45]. Two target genes of miRNAs bra-miRn28-5p and bra-miRn90-5p were PME12 and PME17 (Table S4), and PME activity is important for the virulence of the necrotrophic fungal pathogen Botrytis cinerea by regulate the dicot cell wall pectin content [46]. miRNAs bra-miRn91a-3p and bra-miRn91b-3p shared the same target genes, including CYP79B2, CYP81F3, PUB24, CRK5, ROS1 and LECRK-VII.2 (Table S4). CYP79B2 and CYP81F3 encoding enzymes of indole glucosinolates biosynthesis are coordinately induced in response to P. brassicae [47]. PUB24 acted as negative regulators of PTI in response to several distinct PAMPs [48]. Chen et al. [49] showed that overexpression of CRK5 could trigger hypersensitive response-like cell death in transgenic Arabidopsis plants. Arabidopsis ros1 mutant displayed enhanced resistance to the biotrophic pathogen H. arabidopsidis [50]. LECRK-VII.2 coding a L-type lectin receptor kinase plays a critical role in transfer defense signaling [51]. PINP1, a predicted gene of bra-miR167f-3p (Table S4), regulates RNA silencing pathway in plants [52]. Therefore, we hypothesized a regulatory network of non-heading Chinese cabbage defense response to downy mildew by cell wall modification, metabolite regulation and defense-related gene expression. This study would fills the gaps in miRNAs-mediated downy mildew resistance in nonheading Chinese cabbage and provide some clues to plant miRNA research (Fig. 9).

Conclusions
This research identified miRNAs associated with downy mildew infection in non-heading Chinese cabbage. Target genes and network analysis contributed to understanding the interaction between plant and pathogen. We laid the foundation for further study of miRNA design and transformation to elucidate the function and regulatory mechanisms of plant miRNAs in disease resistance.

Plant materials and treatments
A downy mildew-resistant cultivar (NHCC001) of non-heading Chinese cabbage (B. rapa ssp. chinensis) was chosen as plant material. NHCC001 derived from self-pollinated plants of at least 6 generations in our lab have been re-sequencing by Song at al. [53] (http://nhccdata.njau.edu.cn/). Seedlings were grown in the growth chamber under 25℃ with a 12 h-light/12 h-dark cycle.
Preparation and inoculation of H. parasitica inoculum were performed as described by Sun et al. [54].
Leaves were collected at 0, 6, 12, 24 and 48 hours post inoculation (hpi), frozen in liquid nitrogen and stored at -80℃ until RNA extraction. For each time point, three biological replicates were involved.

RNA extraction
Total RNA was extracted from leaf samples by using a modified CTAB method and treated with DNA erasol to avoid DNA contamination. RNA purity was checked using a NanoPhotometer

Construction and sequencing of Small RNA libraries
A total of 1.5 μg RNA per sample was used as input material for the RNA sample preparations. Fifteen small RNA libraries (3 biological replicates for five time points) were generated using a NEBNext® Ultra™ II small RNA Sample Library Prep Kit for Illumina (NEB, USA) according to the manufacturer's protocol. Because the small RNAs contained 5'phosphate groups and 3'hydroxyl groups, T4 RNA ligase was used to add adaptors at both ends. Libraries were constructed through a series of experiments including reverse transcription, PCR amplification and gel extraction. The quality of each library was assessed on an Agilent Bioanalyzer 2100 system before sequencing, and the small RNA sequencing was performed on an Illunima HiSeq X-ten platform, which generated 50-bp single-end reads. The sequencing work were performed by the Biomarker Technologies Corporation (Beijing, China). The full data sets have been submitted to the Gene Expression Omnibus (GEO) database of NCBI under SRA accession SUB5672046, Bioproject: PRJNA544898.

Analysis of small RNA sequencing data
The sequencing reads were processed to remove adaptors and cleaned by Q30 value. Reads with >20% bases <Q30, and N base >10% were filtered, and sequences shorter than 18 nt or longer than 30 nt were also removed. Silva database, GtRNAdb database, Rfam database and Repbase database were applied to find out ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), and other non-coding RNA as well as repeat sequences by using Bowtie software [55]. The remaining reads were mapped to Brassica rapa genome (v1.5) (http://brassicadb.org/brad/) and analyzed by miRBase (v22) (http://www.mirbase.org) and miRDeep-P [29]. According to the identification method described by Jiang et al. [56], reads that matched to miRBase with no more than 3 mismatches were called as known miRNAs; otherwise, thier were classified into novel miRNAs.

Differential expression analysis of miRNAs
To quantify the abundance of each miRNA, the Transcripts Per Kilobase Million (TPM) value was defined as 'Read count × 1,000,000'/'Mapped Reads'. Differential expression analysis of miRNAs were performed using the DESeq2 R package (1.10.1) with default parameters [31]. DESeq2 provided statistical routines for determining differential expression in digital miRNA expression data using a model based on the negative binomial distribution. The resulting P values were adjusted by using the Benjamini and Hochberg's approach for controlling the false discovery rate. miRNA with |log2(FC)| ≥1.00 and FDR≤0.05 found by DESeq2 were assigned as differentially expressed.

Prediction of potential targets
The software TargetFinder (1.6) [57] were employed to predicted the miRNA targets based on sequence information of known miRNAs and novel miRNAs in Brassica rapa genome (v1.5) (http://brassicadb.org/brad/). The number of mismatches at complementary sites, where there was no gap, between miRNA sequence and potential mRNA target were less than 4 [58]. The miRNA targets were confirmed using psRNATarget online server (http://plantgrn.noble.org/psRNATarget/). The entire set of miRNA target genes was annotated against the NR (ftp://ftp.ncbi.nih.gov/blast/db/), Swiss-Prot

GO and KEGG pathway enrichment analysis
Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GOseq R packages based on Wallenius non-central hyper-geometric distribution.
KEGG [59] is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other highthroughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS [60] software to test the statistical enrichment of differential expression genes in KEGG pathways.

qRT-PCR validation
cDNA was synthesized from 1 µg total RNA with the HiScript ® II Q RT SurperMix for qPCR (+gDNA wiper) reagent Kit (Vazyme Biotche Co., Ltd) and diluted 1/10 by ddH 2 O. qRT-PCR assay was carried out with 2 µL of cDNA and 10 µL of 2×TransStart ® Tip Green qPCR SuperMix (Vazyme Biotche Co., Ltd) using the 7500 Fast Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Actin-7 was chosen as the normalizer and data were analyzed by using 2 -ΔΔCt method [61]. Three technical replicates were ran for each sample. The primers used for qRT-PCR are listed in Table S1

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Funding
The financial of this research was supported by National Natural Science Foundation of China (31471886 and 31272173). The funders had no role in the study design, data analysis and collection, data interpretation and manuscript preparation.

Authors' contributions
performed the experiments. CZS, DH, YW and H-WT analyzed the data. C-ZS, DH, YW and YL wrote and revised the paper. All authors read and approved the final version of manuscript.  Table 2 Classification of small RNAs filtrated.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download. Table S4.xlsx   Table S1.xlsx Table S2.xlsx Table S3.xlsx