Combined small RNA and gene expression analysis revealed roles of miRNAs in maize response to rice black-streaked dwarf virus infection

Maize rough dwarf disease, caused by rice black-streaked dwarf virus (RBSDV), is a devastating disease in maize (Zea mays L.). MicroRNAs (miRNAs) are known to play critical roles in regulation of plant growth, development, and adaptation to abiotic and biotic stresses. To elucidate the roles of miRNAs in the regulation of maize in response to RBSDV, we employed high-throughput sequencing technology to analyze the miRNAome and transcriptome following RBSDV infection. A total of 76 known miRNAs, 226 potential novel miRNAs and 351 target genes were identified. Our dataset showed that the expression patterns of 81 miRNAs changed dramatically in response to RBSDV infection. Transcriptome analysis showed that 453 genes were differentially expressed after RBSDV infection. GO, COG and KEGG analysis results demonstrated that genes involved with photosynthesis and metabolism were significantly enriched. In addition, twelve miRNA-mRNA interaction pairs were identified, and six of them were likely to play significant roles in maize response to RBSDV. This study provided valuable information for understanding the molecular mechanism of maize disease resistance, and could be useful in method development to protect maize against RBSDV.


Results
Small RNA deep sequencing and data analysis. Maize B73 was naturally infected by RBSDV in the field where the Maize Rough Dwarf Disease happened seriously. The control plants were grown in the field and covered by net to prevent the planthoppers, the carrier of the virus. To test whether plants were infected by RBSDV, two pair of primers pS6-604 and pS7-342 were designed according to the genome sequences of segment S6 (GenBank No: HF955010) and segment S7 (GenBank No: HF955011) of RBSDV, respectively. These two pair of primers were used to amplify in the maize individuals with RBSDV symptoms. As a result, the virus genes were amplified in all treatment plants, and could not be detected in control plants ( Supplementary Fig. S1). These results suggested that the phenotype/symptoms were caused by RBSDV. To identify small RNAs from maize, two libraries generated using RBSDV infected plants (TL1 and TL2) and two libraries (CL1 and CL2) generated using the control plants were constructed for high-throughput sequencing. A total of 23,056,821, 20,003,963, 26,678,964 and 20,585,338 raw reads were obtained from CL1, CL2, TL1 and TL2, respectively (Supplementary  Table S1). After removing low quality reads, reads less than 18 nt and reads longer than 29 nt, a total of 9,274,931, 8,351,418,9,552,470 and 13,037,763 clean reads remained from CL1, CL2, TL1 and TL2 libraries, respectively (Supplementary Table S1). These clean reads were used for further analysis. Firstly, clean reads were aligned with maize genome (B73 RefGen_V2, release 5b.60) and Rfam database. Reads annotated as rRNA, snRNA (small nuclear RNAs), snoRNA (nucleolar RNAs), repbase (reads positioned at repeat loci) and tRNA were identified (Supplementary Table S2). The distribution of small RNAs identified from CL and TL libraries is summarized (Fig. 1a). It was also shown that 21-nt small RNAs were the predominant class in maize, followed by 22-nt and 24-nt small RNAs. After infected with RBSDV, the number of 20-, 21-and 22-nt small RNAs in TL libraries increased, while the number of other small RNAs decreased compared with CL libraries. The first nucleotide bias of small RNAs was analyzed. For the small RNAs of 20-22 nt, the canonical length of miRNAs, a strong bias for U of the first nucleotide was observed (Fig. 1b). The small RNA sequencing data has been deposited in NCBI Short Read Archive (SRA) database (BioProject ID: PRJNA438075, Accession number: SRR6829172-SRR6829175).

Identification of known and novel miRNAs in maize.
Comparison the maize small RNAs to the miRBase allowed us to identify 76 known miRNAs, belonging to 26 miRNA families (Supplementary Table S3). Previous studies found that there were twenty miRNAs families were conserved in Arabidopsis, Oryza ostiva and Populus trichocarpa 13,24 . In our dataset, nineteen conserved miRNA families were detected in maize. In addition, seven known but non-conserved miRNA families including MIR408, MIR528, MIR529, MIR827, MIR1432, MIR2118 and MIR2275 were also identified. The normalized expression level of miR166f was 103,108 (TP10M, Numbers of tags per ten million), representing the most abundant miRNAs. The abundance of miR395o-5p, miR2118d, miR395k-3p, miR167g-3p, miR167c-3p, miR408b-3p, miR169r-3p, miR2275a-5p, miR398b-3p and miR398a-3p was low in both sRNA libraries (Supplementary Table S3). According to the criteria as described in previously 25 , a total of 226 potential novel miRNAs were identified. The length of novel miRNAs ranged from 20 to 22 nt, and more than 93% novel miRNAs with the length of 21-22 nt (Supplementary Table S4). The negative folding free energies of these precursors hairpin structures ranged from −89.8 to −16.1 (kcal/mol) with an average of −44.08 kcal/mol, which is similar to the results from other plants. Some of these novel miRNAs were specifically detected in control or treatment libraries. For example, zma-miRn177 was detected only in control libraries, while zma-miRn223 and zma-miRn224 were observed only in treatment libraries.

Target prediction of maize miRNAs and function annotation.
To gain a better understanding of the regulation roles of maize known and novel miRNAs, target genes were predicted using psRNATarget software by comparing miRNA sequences against maize B73 reference genome. A total of 213 targets of 50 known miRNAs and 138 targets of 75 novel miRNA candidates were identified (Supplementary Table S5-S6). Functional annotation of these target genes showed that many defense related genes were regulated by known miRNAs. For example, peroxidase 2 gene (GRMZM2G427815), LRR receptor-like serine/threonine-protein kinase gene (GRMZM2G304745), and a Zea mays rust resistance protein rp3-1 (GRMZM2G045955) were targeted by zma-miR399a-5p, zma-miR390b-5p and zma-miR528b-5p, respectively (Supplementary Table S5). For the targets of novel miRNAs, 37.68% of them were genes with unknown function and 20% of them were transcription factors. Our data showed that many defense related genes were also regulated by novel miRNAs including a plant viral-response family protein gene (GRMZM2G171036), which was regulated by miRn200 (Fig. 2, Supplementary Table S6).

MiRNA expression profiles upon RBSDV infection.
To analyze the expression change of miRNAs in response to RBSDV, the abundance of miRNAs was normalized using numbers of tags per ten million (TP10M), and the relatively expression level of miRNA was calculating by log 2 Ratio (TL/CL). A total of 81 differentially expressed miRNAs (|log 2 | ≥ 1.0) were identified, including 26 known miRNAs and 55 novel miRNAs (Fig. 3). Interestingly, the 26 differential expressed known miRNAs were all up-regulated, and the overall expression levels of all known miRNAs showed up-regulation trend in response to RBSDV (Fig. 4, Supplementary Table S7). Among the 55 differential expressed novel miRNAs, 45 were up-regulated and 10 were down-regulated (Fig. 3). MiRn177 was expressed only in control samples, miRn223 and miRn224 expressed only in virus treated samples. In addition,  our results showed that the expression levels of four novel miRNAs, miRn222, miRn225, miRn226 and miRn146 were significantly increased after infected with RBSDV ( Supplementary Fig. S2, Supplementary Table S4).
Global mRNA expression profiles in maize in response to virus infection. In order to identify the global gene expression alteration upon virus infection, we used the next generation sequencing technology to analyze the transcriptome of maize before and after virus infection. A total of 8.13 Gb data was generated, comprised of more than 40 million reads (Supplementary Table S8). Sequencing randomness analysis was tested for estimating the gene whether or not random distributed in different positions on each genes. The statistical analysis result showed that the sequencing in all samples was in good randomness ( Supplementary Fig. S3a). Saturation analysis showed that the number of genes increased with the total number of tags and reached a plateau after 2.5 million tags ( Supplementary Fig. S3b). Using TopHat alignment, more than 89% of the reads could be successfully mapped to B73 genome, which covers 27,554 genes. The RNA-seq data have been deposited at SRA database (BioProject ID: PRJNA438075, Accession number: SRR6829168-SRR6829171).

Identification and function analysis of differentially expressed genes (DEGs).
To identify differential expressed genes (DEGs) response to virus infection, comparison analysis between control and treated transcriptome libraries was performed. The expression level of genes were normalized by FPKM (expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced). The pearson correlation values between two control (E1 and E3) and two treatment libraries (E2 and E4) were 0.964 and 0.948, respectively ( Supplementary Fig. S3c,d). Under the criterion of P-value ≤ 0.001 and |log 2 | ≥ 1.0, a total of 453 DEGs were found, including 260 up-regulated and 193 down-regulated genes (Fig. 5). The function of these genes were annotated by alignment with Nr and SWISSPROT Database (Supplementary Table S9). Functional annotation indicated that many disease resistance related genes were up-regulated after RBSDV infection. For example, glutathione S-transferase, lipoxygenase, lectin-like receptor protein kinase, O-methyltransferase 8, pathogenesis-related protein 10 and non-specific lipid-transfer protein, etc. (Supplementary Table S9).
In order to explore the functions of these differential expressed genes that are responsive to RBSDV infection, Gene ontology (GO), COG annotation and Pathway enrichment analysis were performed. From our dataset, 168 of 453 differential expressed genes have significant homologies in COG database and were assigned into 25 categories (Supplementary Table S9, Supplementary Fig. S4). Among them, "General function prediction only", "Carbohydrate transport and metabolism", "Replication, recombination and repair", "Amino acid transport and metabolism" and "Energy production and conversion" ranked the top five categories (Supplementary Table S9, Supplementary Fig. S4). To further understand the metabolic and regulatory process for RBSDV-responsive, all upand down-regulated genes were subjected to BGI WEGO program for GO analysis. The detailed summary of GO classification showed that cell, cell part, organelle were the most abundant ones in cell component categories. About molecular function category, the most abundant were binding and catalytic. The last category was biological process, in which cellular process, metabolic process, and response to stimulus were enriched ( Supplementary Fig. S5).
According to KEGG analysis, 77 DEGs were annotated into 65 pathways. Among them, 23 pathways were enriched in response to RBSDV infection including two photosynthesis related pathways, photosynthesis and carbon fixation in photosynthetic organisms (Table 1). In addition, many DEGs were enriched in pathways involved in metabolite or secondary metabolite synthesis, such as starch and sucrose metabolism, pyruvate metabolism,  butanoate metabolism, alanine, beta-Alanine metabolism, glutathione metabolism, aspartate and glutamate metabolism, suggesting significant metabolic changes after RBSDV infection. We found four genes, including a putative coronatine-insensitive protein (GRMZM2G035314, log 2 = 2.16), a respiratory burst oxidase-like protein (GRMZM2G043435, log 2 = 1.07), a heat shock protein HSP82 (GRMZM5G833699, log 2 = 1.78) and one unknown protein (GRMZM2G151519, log 2 = 1.44), were all enriched in plant-pathogen interaction pathway. Interestingly, these genes were all up-regulated in response to virus infection (Supplementary Table S9).

The expression of defense related genes response to RBSDV infection. Previous research has
demonstrated that the expression of some defense response-related genes changed significantly when plants were suffered with biotic or abiotic stresses 7,26,27 . In this study, we investigated the expression changes of defense related genes, and found the expression of 90 defense related genes altered after RBSDV infection, including 53 receptor-like protein kinase genes, six WRKY DNA-binding protein genes, two NBS-LRR family genes, nine pathogenesis-related genes, 13 glutathione S-transferase genes, three peroxidase genes, one heat shock protein  gene, one ferredoxin and two other disease resistance analog genes ( Table 2). Receptor-like protein kinase was an important signal introduction factor involved in plant disease resistance 28 . As shown in Table 2, more than half of these defense related genes encoded diverse receptor-like protein kinases. Among them, the expression of 39 receptor-like protein kinase were up-regulated, while 14 of them were down-regulated upon virus infection. We found that the expression of some of these genes were altered dramatically, for example, the glutathione S-transferase gene (GRMZM2G146913) was up-regulated for about 30 times (log 2 = 4.9274), the peroxidase gene (GRMZM2G313184) was up-regulated for about 24 times (log 2 = 4.6027). These results provided important information for us to understand the mechanism under miRNA regulation of disease resistance in maize.

Quantitative real-time PCR validation.
To validate the deep-sequencing data, we used quantitative real-time PCR (qRT-PCR) to analyze the expression of miRNAs and mRNAs. Ten miRNAs were selected for qRT-PCR analysis including six known miRNAs and four novel miRNAs. Ten genes were also selected for qRT-PCR analysis. These genes included eight genes related with stress response and two genes related with hormone synthesis and metabolism (Fig. 6, Supplementary Table S10). Pearson correlation values between qRT-PCR and RNA-seq with R = 0.824, suggesting that the sequencing data was consistent with the qRT-PCR results.

Combined expression analysis of miRNAs and their targets after virus infection. In recent years,
high-throughput sequencing method has become a powerful technology for global transcriptome and miR-NAome analysis. It has been widely used in many plant species. Here, we simultaneous analyzed the miRNA and gene expression using the same samples before and after RBSDV infection. Through analyzing the relationship between miRNAs and their target genes, twelve miRNA-mRNA pairs were identified, which showed opposing expression patterns response to virus infection (Table 3). In rice, miR156/miR529 and SQUAMOSA PROMOTER BINDING LIKE PROTEIN (SPL) genes constituted a spatiotemporally coordinated gene network which playing an important regulation roles in tiller and panicle branching 29,30 . Plant SPL genes were involved in leaf development, gibberellin response, light signaling, copper homeostasis, response to stresses, and positively regulate inflorescence meristem 29 . We found miR529 were up-regulated in the virus infected samples, and the expression of its three target genes (SPL genes) were all down-regulated after infected with virus in maize (Table 3, Fig. 7). These results were coincided with the significant phenotype changes of virus infected maize, including the abnormal leaf morphology, dwarf, and the abnormality in vegetative and reproductive growth. One target gene of miR528, which was down-regulated after RBSDV infection, is highly homologous with maize rust resistance protein rp3-1, an important defense gene in maize rust resistance caused by Puccinia sorghi 31 . One of the target gene of miR408b-3p is endoglucanase, which catalyzes the hydrolysis of cellulose. Our results showed that miR408b-3p was up-regulated significantly, while its target gene was down-regulated upon virus infection (Table 2, Fig. 7). In addition, the expression trend  of miR399d-5p, miR398b-5p and miR156k showed a negative correlation with their target genes. Three novel miRNA-mRNA pairs were identified, of which GRMZM2G076468 encodes a cyclin-dependent protein kinase (Cyclin-P4-1), is the target of zma-miRn53. GRMZM2G160917 is a SPL gene and regulated by zma-miRn138. Zma-miRn194 targeted GRMZM2G484653, a gene annotated as unknown function ( Table 2, Fig. 7).
Accumulated evidence indicated that many plant endogenous miRNAs were responsive to pathogenic fungus and virus infection and played important roles in plant disease resistance, such as miR398, miR160 16 , miR482, miR2118 17 , and miR472 18 . After infection with Magnaporthe Oryzae, the expression of rice miR398 was induced both in susceptible and resistant lines, while miR160 were only induced in resistant lines, and overexpression of miR160a or miR398b can enhance rice resistance to the disease 16 . Here, we found that the accumulation of miR398 were significantly induced to higher levels upon RBSDV infection, and miR160a decreased upon RBSDV infection (Supplementary Table S3, Fig. 7). The expression trend of miR160 and miR398 in susceptible maize variety B73 was consistent with that in susceptible lines of rice after infected with pathogen. These results suggested miR160 and miR398 might played great roles in maize disease resistance. In addition, miR482, miR2118 and miR472 contribute to plant immunity through negative regulation of R gene 17,18 . However, the expression level of miR2118 was very low, while miR482 and miR472 were not detected in our dataset.
The expression of genes that involved in dwarf symptoms. After infected with RBSDV, the growth and development of maize exhibited severe abnormalities, such as dwarf, dark green leaf, and failure of completing the reproductive growth in most cases 2,32 . The height of plants is shown to correlate with the composition of the cell wall, which are associated with metabolism and biosynthesis of cellulose, lignin, hemicellulose and pectin 7,33,34 . Cellulose is a major class of polysaccharide, which is the main ingredients of plant cell wall, and was closely relative to plant defense 35 . We identified a variety of gene families involved in cell wall synthesis and degradation, and the expression of these genes altered when the maize infected by RBSDV (Table 4). These genes included 13 cellulose synthase, seven endoglucanase, two glycoside transferase, one glycosyltransferase, five pectin related genes, three glucosidase, galacturonase, one putative mixed-linked glucan synthase and one glycine-rich cell wall structural protein gene (Table 4). Cellulose synthase is an important enzymes important in cellulose synthesis system. In the present study, nine cellulose synthase genes were up-regulated and three cellulose synthase genes were down-regulated upon virus infection. Endoglucanase is a specific enzyme that catalyzes the hydrolysis of cellulose, and the expression of seven endoglucanase encoding genes was altered, five were up-regulated, and two were down-regulated (Table 4). Beta-glucosidase catalyzes the hydrolysis of glycosidic bonds, and a variety of glycosidic conjugates of hormones and defense compounds can be hydrolyzed by beta-glucosidases. Here, we found three beta-glucosidases were up-regulated upon the RBSDV infection. The gene encoding a beta-glucosidases 18 gene was up-regulated for almost thirty times. The up-regulated expression of beta-glucosidases might lead to the degradation of defense compounds, and then resulting in the collapse of maize defense system. Plant height is regulated by hormones, such as gibberellins (GAs), auxin (IAA) and brassinosteroid [36][37][38] . In maize, genes that have large effects on plant height have been well characterized, and most of them were involved in hormone synthesis, transport, and signaling, for example, brachytic2 38 , dwarf3 39 and nana plant1 40 . In this study, a total of 17 GA biosynthesis and signaling genes were up-or down-regulated upon RBSDV infection, including two DELLA protein genes, four gibberellin 2-beta-dioxygenase genes, seven gibberellin receptor GID1 genes, two gibberellin 20 oxidase (GA20ox) genes, one chitin-inducible gibberellin-responsive gene and one putative response to gibberellin stimulus genes (Table 5). Auxin plays essential roles in regulating plant growth and development, and also regarded as a negative regulator for plant disease resistance 41,42 . In response to RBSDV infection, the expression of many auxin synthesis and transport related genes was alerted. Among them, seven and four auxin responsive factor (ARF) genes were up-and down-regulated upon RBSDV infection, respectively (Table 5). For example, GRMZM2G390641, encoding ARF21 gene, was regulated by zma-miR160d-5p and miRn91 (Supplementary Table S5   Auxin polar transport is essential for the formation and maintenance of local auxin gradients of plant 43,44 . Auxin efflux carriers PIN family genes played important roles in auxin polar transport. Loss of function of PIN genes severely affected organ initiation. For example, the auxin transport-defective mutants br2 and sem1 showed dwarf phenotype and vasculature defects 43 . Here, the expression level of seven auxin efflux carrier genes were decreased after virus infection (Table 5). It is possible that the decreasing expression of auxin efflux carrier protein genes could be a major reason that caused the dwarf phenotype after maize infected by RBSDV.
In conclusion, we identified 302 miRNAs and 351 potential target genes in maize. The expression patterns of 81 miRNAs differed dramatically upon RBSDV infection. Combined small RNA and gene expression analysis identified 12 miRNA-mRNA pairs with opposite expression patterns response to virus infection, and six of them are likely to play significant roles in the formation of maize disease symptoms. This study provided insight into the roles of miRNAs in response to RBSDV, and could help to develop novel strategy for crops against virus infection.

Materials and Methods
Zea mays B73 was planted in the field where the RBSDV disease happened seriously almost every year. As the control, the plants were covered with a net to prevent the planthoppers. Leaves of one-month-old maize seedlings with rough dwarf disease symptoms were collected. Total RNA were prepared separated from each individual sample using RNAiso Plus reagent (Takara, Dalian, China), following by RNase-free DNase treatment (Takara, Dalian, China). RNA concentration was quantified by Eppendorf BioPhotometer plus UV-Visible Spectrophotometer. The cDNA were synthesis using One Step PrimeScript miRNA cDNA Synthesis Kit (Takara, Dalian, China) according the manufacturer's instructions. According to the sequences of RBSDV segment S6 and S7 (HF955010, HF955011), two pair of primers pS6-604 (5′-CCTAGTTCTCCGCAAGCC-3′, 5′-CAGGGACAGTTCCAATCATAAA-3′) and pS7-342 (5′-TCAGCAAAAGGTAAAGGAAGG -3′, 5′-GCTCCTACTGAGTTGCCTGTC-3′) were designed. Samples were collected according the method in previous studies 7,45 . Every ten virus infected seedlings which can be detected by both the primer pS6-604 and pS7-342 were harvested as one sample, and three more samples were prepared by the same method as replicates.
Construction and sequencing of small RNA libraries. Small RNA libraries were constructed as described in the previous studies 15,46 . Briefly, low molecular weight RNAs (10 nt -40 nt) were isolated from the total RNA by electrophoresis using 15% TBE-urea denaturing polyacrylamide gel. Then, the 5′ and 3′ adaptors were added and reverse transcription was performed to synthesize cDNA. And cDNA libraries were sequenced using Illumina HiSeqTM 2000. The sequencing was accomplished by BGI small RNA pipeline (BGI, Shenzhen, China). After sequencing, clean reads were generated by removing the adapter sequences and low quantity reads (reads having 'N' , <18 nt, and >29 nt). Then clean reads were used to align with maize B73 RefGen_V2 genome (http://archive.maizesequence.org/index.html), GenBank and Rfam database, and miRbase (http://www.mirbase. org/). The detail processes to identify known and novel miRNAs were according to the method described in previous studies 47 . Transcriptome sequencing and bioinformatics analysis. Transcriptome libraries were constructed using Illumina sample preparation kits. Briefly, poly A mRNAs were isolated and cut to short fragments. The short mRNA fragments were then used to synthesize the first strand cDNA using random hexamers primers. Then, dNTPs, RNase H, buffer and DNA polymerase I were added to synthesize the second strand cDNA. cDNAs were further purified using QiaQuick PCR kit. Then, polyA tails and adaptors were added and the DNA fragments with suitable size were recovered from gel. Finally, the cDNA were amplified by PCR, followed by sequencing using Illumina HiSeq ™ 2000.
After sequencing, raw data was filtered to generate clean reads by removing adaptor sequences, reads containing multiple N and lower quality reads. Then, the clean reads were used to compare with maize genome sequences (B73 RefGen_V2, release 5b.60) using SOAPaligner/SOAP2 with the parameters that mismatch ≤2 48 . The gene expression level is calculated by using FPKM method 49 . Differential expression analysis of two samples was performed using rigorous algorithm method with P-value ≤ 0.001 and the absolute value of log2Ratio ≥1. Gene function analysis was carried out by BLASTx searches against the UniProt database and the Swiss-Prot protein database (http://www.expasy.ch/ sprot). Gene Ontology (GO) annotation analysis was based on WEGO (http:// wego.genomics.org.cn/cgi-bin/wego/index.pl). Cluster of Orthologous Groups (COG) classification analysis was based on the database (http://www.ncbi.nlm.nih.gov/COG/). Pathway-based analysis was performed according to Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG) database (http://www.genome.jp/kegg/).

qRT-PCR Validation.
For qRT-PCR validation of miRNAs, the Mir-X miRNA qRT-PCR SYBR Kit (Clontech Laboratories, Inc) were used following the manufacturer's instructions. For all miRNAs and genes, the qRT-PCR was performed in ABI7500 (Applied Biosystems). Primers used for qRT-PCR were listed in Supplementary Table S11. All reactions were performed in biological triplicates. For qRT-PCR analysis of miRNAs and mRNAs, U6 RNA and ubiquitin were used as the internal control, respectively. The relative expression of all mRNAs and miRNAs were calculated using 2 −ΔΔct method.