http://orcid.org/0000-0002-1584-7170
Figures
Abstract
Plant microRNAs (miRNAs) are a class of non-coding RNAs that play important roles in plant development, defense, and symptom development. Here, 547 known miRNAs representing 129 miRNA families, and 282 potential novel miRNAs were identified in Beta macrocarpa using small RNA deep sequencing. A phylogenetic analysis was performed, and 8 Beta lineage-specific miRNAs were identified. Through a differential expression analysis, miRNAs associated with Beet necrotic yellow vein virus (BNYVV) infection were identified and confirmed using a microarray analysis and stem-loop RT-qPCR. In total, 103 known miRNAs representing 38 miRNA families, and 45 potential novel miRNAs were differentially regulated, with at least a two-fold change, in BNYVV-infected plants compared with that of the mock-inoculated control. Targets of these differentially expressed miRNAs were also predicted by degradome sequencing. These differentially expressed miRNAs were involved in hormone biosynthesis and signal transduction pathways, and enhanced axillary bud development and plant defenses. This work is the first to describe miRNAs of the plant genus Beta and may offer a reference for miRNA research in other species in the genus. It provides valuable information on the pathogenicity mechanisms of BNYVV.
Citation: Liu J-Y, Fan H-Y, Wang Y, Zhang Y-L, Li D-W, Yu J-L, et al. (2017) Characterization of microRNAs of Beta macrocarpa and their responses to Beet necrotic yellow vein virus infection. PLoS ONE 12(10): e0186500. https://doi.org/10.1371/journal.pone.0186500
Editor: Sebastien Pfeffer, Institut de Biologie Moleculaire et Cellulaire, FRANCE
Received: April 13, 2017; Accepted: October 3, 2017; Published: October 16, 2017
Copyright: © 2017 Liu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The raw data of deep sequencing and microarray analysis have been deposited in NCBI database. The accession codes of them are GSE101436 and GSE102330.
Funding: This work was supported by the China Agriculture Research System (CARS-170304), www.moa.gov.cn, to CGH, which played a role in the decision to publish and preparation of the manuscript. It was also supported by the National Natural Science Foundation of China (31401708), http://www.nsfc.gov.cn/, to YW, which also participated in the decision to publish. The National Natural Science Foundation of China (31470253), http://www.nsfc.gov.cn/, to YLZ, also played a role in the decision to publish.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Plant microRNAs (miRNAs) are a class of 20–24 nt endogenous small non-coding RNAs. Their coding genes possess their own transcriptional units [1] and are mostly located in intergenic regions, introns or inverse repeats of coding sequences [2]. These miRNA genes (MIR) are transcribed into primary miRNAs with a secondary stem-loop structure in the partial sequence by RNA polymerase II [3, 4]. Primary miRNAs are processed into mature miRNAs through two or more cleavage events, transported to the cytoplasm, methylated and despiralized by helicase. The mature miRNA guide strand enters into the RNA-induced silencing complex and regulates gene expression by target cleavage, translational inhibition and DNA methylation. The miRNA passenger strand (miRNA*) is degraded by an unclear mechanism [5, 6], and some miRNA* may also function in plant defense and development [7, 8].
There are conserved miRNAs and lineage-specific miRNAs in all plant species [9]. miRNA sequence abundance and conservation are correlated [10], and miRNAs are inherited through vertical descent [9]. Some conserved miRNA families, for example miR156, miR160, miR166 and miR172, which are generally highly expressed and ubiquitous across all terrestrial species [10, 11], often have conserved targets that play specific functions in plant development, and are probably under stringent selection pressures [12]. Lineage-specific miRNAs are diverse and common in individual species, even in Arabidopsis thaliana and Arabidopsis lyrata, and are subjected to less selection during evolution within the genus [12, 13].
Genome-wide analyses of many plants infected by viruses have proven that plant viruses often disturb host miRNA pathways, and the aberrant expression of some miRNAs could affect plant defenses or morphological development. For example, Rice stripe virus (RSV) blocks the defense response by significantly down-regulating the miR1870-5p- and miR1423-5p-mediated disease resistance pathways in infected rice plants [14]. Similarly, the induction of miR319 in rice by Rice ragged stunt virus suppresses jasmonic acid (JA)-mediated defenses to facilitate viral infection and symptom development [15]. miR159 plays an important role in eliciting the abnormal phenotype in Arabidopsis infected by Cucumber mosaic virus (strain FNY) [16]; Turnip mosaic virus-infected Arabidopsis misregulates miR167, and the down-regulated miR167 and up-regulated target Auxin Response Factor 8 are major causes for the phenotypic aberrations [17]. Thus, disturbance of the miRNA pathway is closely related to the pathogenesis of plant viruses.
Beet necrotic yellow vein virus (BNYVV) [18], a member of the Benyvirus genus, severely impacts sugar beet production by causing beet rhizomania disease. BNYVV is a multiparticle virus with four or five positive-sense, single-stranded RNAs (RNA1–5) which are individually packaged into rod-shaped virion [19]. The larger RNA1 and RNA2 genomes encode the housekeeping genes of the virus, while the smaller RNAs3, 4, and 5 are involved in pathogenicity and vector transmission during BNYVV systemic infection of beet [20–23]. The RNA3-p25 gene encodes virulence and avirulence factors [21, 24], and amino acid changes in p25 can generate resistance-breaking variants [24]. RNA3-p25 can target the sugar beet 26S proteasome involved in the induction of hypersensitive resistance responses through interactions with an F-box protein [25]. Additionally, other P25-interacting sugar beet proteins represent putative viral targets or components of plant resistance, but their specific functions are not clear [26]. RNA3-p25 can also induce hormonal changes and a root branching phenotype in transgenic Arabidopsis thaliana [27]. In addition, a transcriptome analysis of Beta macrocarpa was performed to identify differentially expressed transcripts in response to BNYVV infection. Gene ontology (GO) analysis showed that these differentially expressed genes were mainly enriched in biotic stimuli and primary metabolic processes [28]. However, to date, miRNAs expression profiles in response to BNYVV infection remains unexplored. In addition, sugar beets, one of the main sugar-yielding crops, are of significant economic value, and miRNA-based biotechnology plays an important role in plant improvement [5]. Although whole-genome information of Beta vulgaris was available in 2014, miRNAs of the plant species in the genus Beta remains largely unknown.
In this study, small RNA deep sequencing was performed to investigate miRNAs from leaves of BNYVV-infected (VL) and virus-free (L) plants. The characteristics of miRNAs were systematically studied through a phylogenetic analysis, and the differential expressions of miRNAs in response to BNYVV infection were investigated and further validated by microarray analysis and stem-loop RT-qPCR. miRNAs’ targets were also predicted by degradome sequencing. To our knowledge, this is the first report that defines B. macrocarpa miRNAs, which would expand our knowledge of miRNAs in the genus Beta, and aid in the discovery of miRNAs related to the pathogenic mechanisms underlying beet–BNYVV interactions.
Results and discussion
High-throughput small RNA sequencing of BNYVV-infected and virus-free leaves of B. macrocarpa
To investigate the B. macrocarpa miRNAs involved in BNYVV infection, two small RNA libraries were constructed from leaves of virus-free (L) and BNYVV-infected (VL) plants. The raw sequence data have been deposited in NCBI database with the accession number GSE101436. Small RNA sequencing generated 17,249,248 and 13,628,866 raw reads from the L and VL libraries, respectively (Table 1). After removal of low-quality and corrupted adaptor sequences (reads < 18 nt), 14,556,099 and 11,458,601 clean reads were obtained from the L and VL libraries, respectively (Table 1). The size distribution of sequences showed that the majority of small RNAs in our libraries belonged to the 24-nt size class (Fig 1A and 1B). In the small RNA libraries from BNYVV-infected plants, 15.36% of the high-quality reads and 5.09% of the unique reads were well matched to the BNYVV genome, whereas negligible numbers (less than 50 reads) were obtained from the healthy control (Table 1). Our sequence data are in agreement with those of many plant species, such as Nicotiana benthamiana [29], Medicago truncatula [30], Punica granatum [31] and Arabidopsis [32] in which the 24-nt size class is the most abundant of the small RNAs. Thus, these data confirmed the high quality of the small RNA libraries, indicating that they were suitable for further comparative analyses.
The size distribution of total (A) or unique (B) beet-derived small RNAs and total (C) or unique (D) miRNAs in two libraries. L, small RNA library from leaves; VL, from BNYVV-infected leaves; nt, nucleotide.
We then identified known miRNAs by mapping these unique small RNAs to miRBase v21.0 (http://microrna.sanger.ac.uk), allowing up to two mismatches during the alignment. Others were mapped to transcripts of B. macrocarpa [27] or reference genomes of B. vulgaris in NCBI to identify their precursors, and then each miRNA whose precursor met the criterion [minimal free energy index (MFEI) > 0.9] was defined as a potential novel miRNA. Finally, 547 miRNAs, belonging to 129 known miRNA families, and 282 potential novel miRNAs, with three or more reads in the L or VL library were identified (S1 Table). Detection of the corresponding miRNA* is an important factor in verifying the existence of miRNA. miRNA-3p and -5p were detected simultaneously in 16 of 129 known miRNA families (S2 Table) and 43 of 282 potential novel miRNAs (S3 Table). However, due to the incomplete genomic information for B. macrocarpa, pre-miRNAs of 97 of the 129 known miRNA families were not identified, including conserved miRNA 168, and miR390 (S1 Table). Thus, all of the known miRNAs identified in this study, with or without proper pre-miRNAs, were considered in the following analysis. The statistical data on the length distribution of mature miRNAs revealed that 21-nt and 22-nt miRNAs were the most abundant miRNAs, and the species of 21-nt and 24-nt miRNAs were the most abundant unique miRNAs (Fig 1C and 1D). We also determined the frequency of the first base of mature miRNAs and found that the 20-nt, 21-nt, and 22-nt miRNAs preferentially started with ‘U’ (73.14%, 83.01% and 94.09%, respectively), while 24-nt miRNAs preferred ‘A’ at the first base (Fig 2).
Phylogenetic analysis of the microRNAome of B. macrocarpa
To gain a deeper understanding of the B. macrocarpa microRNAome, and to offer a reference for the microRNAome research of other species of the genus Beta, or even of Caryophyllales, the B. macrocarpa miRNAs’ distribution within the phylogenetic perspective of plant miRNAs proposed by Taylor et al. [9] was investigated. Through sequence analyses, we identified 28 out of 34 conserved miRNA families in Eudicotyledons, including miR482, miR535, miR536, miR828, miR3627 and miR3630 (Fig 3). They were confirmed by the detection of their corresponding miRNA*s, and a microarray analysis. In total, 13 out of 28 miRNA families were further confirmed by these analyses (S2 Table, S5 Table in bold font). Additionally, 10 of 15 other miRNA families were tested again by microarray analysis (S5 Table in bold font). miR2111, miR397, miR399, miR428 and miR828 were not tested again by microarray analysis because of their relatively low levels of expression, and the precursors of the first three miRNAs met the criteria (MFEI > 1) (S1 Table). No precursors were found for the latter two miRNAs. In healthy B. macrocarpa, the miR159 family with 22 members presented the highest expression abundance with 71,771 reads, followed by miR393. miR535 was expressed at a relatively high level, with more than 10,000 reads compared with its low expression in other species.
The phylogenetic perspective was proposed by Taylor et al. [9]. The graphic mode has been used according to the model of Yin et al. [33]. Black underlined numbers indicate miRNA families in beet; Black italic numbers indicate missing miRNA families. Black * represents the existence of the miRNA was confirmed by detection of its corresponding miRNA*; and red * represents the expression of miRNA was confirmed by a microarray analysis. miR1515 is highlighted in blue, suggesting that it may have been generated before the separation of Citrus and Beta.
There are six missing miRNAs including miR477, miR530, miR2950, miR827, miR4376, and miR4414 (Fig 3). However, there were also no these missing miRNAs in Silene latifolia belonging to the Caryophyllaceae family of Caryophyllales, [10]. We cannot conclude that this is a common characteristic of plant miRNAs in Caryophyllales because of the limited number of samples.
We also found 93 reads of miR1515 (with only one mismatch), a Citrus sinensis species-specific miRNA [9], in the B. macrocarpa microRNAome by small RNA sequencing and a microarray analysis (S5 Table). Thus, miR1515 may have originated before the separation of Citrus and Beta (Fig 3).
To gain more evolutionary insights into the B. macrocarpa microRNAome, a phylogenetic analysis of precursors of miR160 (MIR160) with 29 species (Fig 4) used in Yin et al. [33] and a cross-species conservative analysis of known miRNAs were performed. To construct the MIR160 phylogenetic tree, ath-miR160a-5p (reads = 871) was chosen as the study object and renamed as bma-miR160a. bma-miR160a was the most abundant of three miR160s found in our small RNA sequencing, and, based on the secondary structure of its precursor, it had the canonical stem-loop structure of a processed mature miRNA (MFEI = 1.10) (S1 File). In addition, bma-miR160a* (reads = 74), which matched its precursor, was also identified. The phylogenetic tree showed that bma-MIR160a formed a subgroup with Glycine max, M. truncatula, A. lyrata, and A. thaliana (Fig 5), and the relationship between bma-MIR160a and gma-, mtr-, aly- or ath-MIR160 was close.
The phylogenetic tree was built using the common tree tool in NCBI.(https://www.ncbi.nlm.nih.gov/Taxonomy/CommonTree/). Species are presented in order using the three-letter codes to the right of the full species name.
The phylogenetic tree was constructed using precursors of miRNA160 from diverse plant species obtained from miRBase v21.0 by MEGA 5.0.
The statistical data of the cross-species conservative analysis of known miRNAs was consistent. When we aligned unique sequences to miRBase v21.0 to identify known miRNAs, these unique sequences sometimes matched precursors of multiple miRNAs of multiple species with no more than two mismatches, and these were termed as the cross-species conservative miRNAs of the corresponding B. macrocarpa miRNAs. Thus, we analyzed the similarities between the miRNAs of B. macrocarpa and those of other species. The top three species having the greatest number of similar miRNAs were G. max, Populus trichocarpa and Malus domestica (227, 164 and 159, respectively) belonging to the fabids. Additionally, there were 127 miRNAs of M. truncatula, 104 of A. thaliana and 97 of A. lyrata that were similar to those of B. macrocarpa (Fig 6A).
Numbers of cross-species conservative miRNAs with (A) or without (B) conserved miRNAs were counted. Conserved miRNAs include miR156, miR159, miR160, miR162, miR164, miR166, miR167, miR168, miR169, miR171, miR172, miR319, miR390, miR393, miR394, miR395, miR396, miR397, miR398, miR399, miR403 (in Eudicotyledons), and miR408.
Most of these miRNAs were conserved miRNAs [11] that exist widely in the Eudicotyledons. Additionally, the greatest numbers of non-conserved miRNAs were found in M. truncatula, G. max and A. thaliana (45, 24 and 16, respectively) (Fig 6B). In conclusion, the relationship between the B. macrocarpa microRNAome and M. truncatula, G. max or A. thaliana, belonging to the rosids, may be close. By contrast, the relationship between the B. macrocarpa microRNAome and Nicotiana tabacum, Solanum lycopersicum or Solanum tuberosum, belonging to the Solanaceae, may be distant because they have few miRNAs that are similar to those of B. macrocarpa (2, 0 and 7, respectively) (Fig 6B), and there are no miRNAs enriched in the Solanaceae as in B. macrocarpa, such as miR1919, miR4376 and miR5300 [10].
In summary, we have uncovered useful information for miRNA research in other species of the genus Beta and for plant evolutionary analyses. There were also some miRNAs of Poaceae plants, most with low expression levels, that were similar to B. macrocarpa miRNAs, but the relationship between Monocotyledons and Eudicotyledons is complex.
The microRNAome of B. macrocarpa has species specificity. The lineage-specific miRNAs in B. macrocarpa were investigated by re-evaluating potential novel miRNAs (reads > 20 in L/VL library) with -3p and -5p according to the stringent criteria reported by Taylor et al. [9]. To validate the lineage-specific miRNAs, we aligned the precursor sequences to the public database in NCBI, and we found that no miRNA gene was supported by other plant species data sets, but most of them were found in B. vulgaris (Table 2). Finally, 8 possible lineage-specific miRNAs were identified, which could have evolved within the genus Beta (S2–S5 Files).
The differential expression analysis of miRNAs responsive to BNYVV infection was confirmed by a microarray analysis and stem-loop RT-qPCR
Genome-wide profiling of miRNAs revealed that 103 known miRNAs (with reads > 20 in L or VL library), representing 36 miRNAs families, and 45 novel miRNAs were differentially regulated by at least a two-fold change in BNYVV-infected B. macrocarpa compared with the mock-inoculated control (|log2(VL/L)| > 1.0). The expression levels of 20 known and 20 novel miRNAs were inhibited during BNYVV infection (S4 Table; S6 File).
To determine the reliability of small RNA sequencing, a microarray analysis was performed with 158 probes of conserved miRNAs from the phylogenetic analysis (80), other known miRNAs (18), and potential novel miRNAs (60) with reads > 20 in L/VL library (S5 Table). The raw sequence data have been deposited in NCBI database with the accession number GSE102330. The expression of 94 miRNAs (hybrid signal > 500 in L or VL microarray library) was further analyzed. Results showed that the fold change of expression obtained by microarray analysis was not completely consistent with small RNA sequencing results, mainly due to the different principles and sensitivity of the two technologies, but the trend of 39 out of 55 significantly differential expressed miRNAs was similar (S5 Table highlighted in yellow). And the trend of the other 16 miRNAs was inconsistent (S5 Table highlighted in red) between microarray analysis and small RNA sequencing. We think there were two possible reasons for such phenomenon. One is that microarray analysis hardly accurately distinguishes the miRNA sequences at the ends of the bases, such as miR159_R+1_1ss21AT, miR159a_R-1_1ss1TG and miR156t, and the other is that two kinds of strong fluorescent signals might mask the difference between each other, such as miR166a_1ss21CT, miR166a and miR535d_1ss7CT.
Meanwhile, the stem-loop RT-qPCR was performed to validate the small RNA deep sequencing results. The expression patterns of eight known miRNAs (miR156, miR166, miR168, miR171, miR172, miR393, miR535 and miR1515) and three novel miRNAs (miRn31-5p, miR111-3p and miRn187-5p) were tested. A high correlation (R2 = 0.85) was found between stem-loop RT-qPCR and small RNA sequencing (Fig 7A and 7B).Thus, the alterations in miRNA expression detected by small RNA sequencing could reflect the actual miRNA expression changes between BNYVV-infected and mock-inoculated control plants. miRNAs were selectively regulated by BNYVV infection. miR169, miR172 and miR319 were inhibited, and miR156, miR160, miR390, miR162, miR168 and miR393 were induced (Table 3).
(A) Expression patterns of selected B. macrocarpa miRNAs in response to BNYVV infection by stem-loop RT-qPCR. (B) Correlation of the differential expression ratio of selected miRNAs measured by stem-loop RT-qPCR and small RNA-Seq. (C) Expression patterns of the targets of selected miRNAs in response to BNYVV infection by RT-qPCR. PP2A was used as an internal control.
Analysis of the predicted target genes of differentially expressed miRNAs.
Degradome sequencing [34] was performed to investigate the targets of miRNAs. One library was constructed from leaves of BNYVV-infected (VL) plants. Through the degradome sequencing approach, we obtained about 7,678,479 raw reads, of which 37,371 reads were perfectly matched to the beet genome.
After further analysis, 32 target genes for 15 known miRNA families and 4 target genes for 4 novel miRNA were identified (Table 4). The sequences of candidate targets and the results of t-plot were shown in S7 and S8 Files. The cleaved targets have been categorized into five classes (Categories 0, 1, 2, 3 and 4) according to the abundance of degradome tags at the target sites [35]. In our study, we listed Categories 0–3 targets, and omitted Category 4 targets which have only one raw read at the predicted cleavage position. Our results substantiate that conserved plant miRNAs regulate conserved targets at identical target sites in species they exist in [12]. This could be manifested by the miR156’s squamosa promoter-binding-like protein, miR162’s endoribonuclease Dicer homolog 1 and miR393’s transport inhibitor response 1.
To investigate the functions of these miRNAs in response to BNYVV infection, the differential expression patterns of the target genes of miR156, miR166, miR169, miR171, miR172 and miR393 were also examined by RT-qPCR (Fig 7C). And the results showed that the differential expression of these target genes appeared to be negatively related to their respective miRNAs.
In addition, miRNAs regulate gene expression, not only by cleaving targets, but also by inhibiting translation and DNA methylation [6]. Tombusviruses infections enhance the accumulation of miR168, a regulator of ARGONAUTE1 (AGO1) mRNA, and, in parallel, induce the expression of AGO1 mRNA, which has an inhibited translational capacity, and the infection decreases the AGO1 content, resulting in disturbance of its anti-viral function [36].
Roles of the miRNA targets’ gene-mediated expressions in plant defenses.
BNYVV infections also disturbed B. macrocarpa’s defense system. The RNA silencing pathway is an efficient way to protect plants from viral infection. miR162-controlled endoribonuclease Dicer homolog 1, miR168-controlled AGO1 and miR403-controlled AGO2 are positive factors participating in the RNA silencing pathway [36]. In BNYVV infected B. macrocarpa, miR162 and miR168 were induced (VL/L > 1.8), thus impairing the plant’s antiviral capability and enhancing viral infection and proliferation (Table 3).
Interactions of miRNA-target-mediated gene expression levels responsible for plant morphological changes.
BNYVV infections lead to chlorosis, dwarfism and enhanced axillary bud development [28]. miR395 mediated regulation of sulfate accumulation and allocation in Arabidopsis thaliana. Overexpression of miR395 inhibited the process of sulfur assimilation, and led plants display sulfur deficiency symptoms. In BNYVV-infected B. macrocarpa, miR395 was significantly up-regulated, and might play a role in the formation of chlorosis and dwarfism symptoms [37]. The abnormal development of axillary buds may result from the increased abundance of miR156. Up-regulated miR156 can disturb gibberellin levels [38], and exhibit enhanced branching from axillary buds resulting in a bushy appearance [39–42]. Additionally, miR156 was induced by the BNYVV infection and may be responsible for abnormal axillary bud development. miR156 is a switch between vegetative and reproductive growth in plants through its regulation of the ‘miR156-squamosa promoter-binding-like protein–miR172- AP2-like factor’ pathway. The function of the miR172’s target is to inhibit flowering. The up-regulated miR156 (VL/L > 4) and the corresponding down-regulated miR172 (VL/L < 0.6) may retain the vegetative state by delaying the flowering time [43, 44]. Up-regulated miR169 plays a key role in stress-induced early flowering in A. thaliana [45], and in BNYVV-infected B. macrocarpa, miR169 (VL/L < 0.3) was inhibited, delaying flowering, which correlated with the miR156 regulatory pathway.
BNYVV infection altered miRNAs participating in plant hormone synthesis and signal transduction pathways.
The plant hormone signaling pathway plays an important role in plant growth and defense. Auxin is related to plant height, root development and apical dominance. miR160, miR390 and miR393 control plant growth by regulating the auxin signaling pathway [46, 47]. miR390 and miR393 were up-regulated (VL/L > 2.2) by BNYVV infection, and the target of miR393, transport inhibitor response 1, was inhibited. All of these may function in the formation of dwarf symptoms. The JA signaling pathway functions in plant disease resistance and insect resistance. miR319-controlled TEOSINTE BRANCHED/CYCLOIDEA/PCF transcription factors control the biosynthesis of the hormone JA by affecting the expression levels of JA biosynthetic genes and is proposed to coordinate two sequential processes in leaf development: leaf growth, which they negatively regulate, and leaf senescence, which they positively regulate [15, 48]. In BNVV-infected B. macrocarpa, miR319 was inhibited (VL/L < 0.4), which might enhance plant defense by increase the JA content, and might also take part in leaf area reduction.
Conclusions
miRNA expression profiles in B. macrocarpa during viral infection were investigated by small RNA sequencing, and were further validated by microarray analysis and stem-loop RT-qPCR. In B. macrocarpa, 547 known miRNAs, representing 129 miRNA families, and 282 potential novel miRNAs were identified, and a phylogenetic analysis showed 8 lineage-specific miRNAs in the genus Beta. The differential expressions of 103 known miRNAs, representing 38 miRNA families, and 45 potential novel miRNAs had at least a two-fold change in response to BNYVV infection, and most of the miRNAs were involved in the auxin signal pathway, jasmonate biosynthesis, and plant defense. Thus, the BNYVV infection disturbed the miRNA pathway of B. macrocarpa, which may be favorable for viral symptom development. Our results revealed miRNAs involved in the BNYVV infection, which would be helpful for further study of the molecular mechanisms underlying BNYVV-plant interactions.
Materials and methods
Plants, viral inoculations and detection
B. macrocarpa plants were grown in a controlled-environment chamber at 24 ± 1°C with 16 h of illumination and 8 h darkness per day. BNYVV (RNAs 1+2+3+4+5) was a mixture of total RNAs (2 μg) from ‘BN34’ (RNAs 1+2+3+4) preserved in our laboratory [49], and in vitro transcript of RNA5 (1 μg). BNYVV RNAs were then mixed with equal volume of inoculation buffer (50 mM glycine, 30 mM K2HPO4, 1% bentonite, and 1% celite, pH 9.2), and rubbed onto one leaf of B. macrocarpa. Inoculation buffer without the addition of RNAs served as the control. Three leaves per plant were inoculated. At 15 days post inoculation (dpi), western blot analysis of the BNYVV coat protein (CP) was performed using symptomatic leaves of BNYVV-inoculated B. macrocarpa or the healthy control plants (S1 Fig). Rabbit polyclonal antibodies against BNYVV CP had been preserved in our laboratory [28].
Total RNA extraction
Total RNAs were extracted from systemic leaves of virus-free and BNYVV-infected B. macrocarpa at 15 dpi using according to previously described methods [28]. RNA integrity and size distribution were examined using a Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA) and agarose gel electrophoresis (1%). For each group, the RNA pool was prepared by mixing RNA samples (12 μg per sample) from five individual plants.
Construction and sequencing of small RNA libraries
Small RNA fragments (10–40 nt) were isolated from total RNA sample (2 μg), and two respective small RNA libraries were constructed using IlluminaTruSeqTM Small RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. 5he performdaptors was added to the short RNA fragments of each sample followed by RT-PCR. The products were purified and assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies). A clustering of two samples was performed using a cBot Cluster Generation System using TruSeqSE Cluster Kit (Illumina) according to the manufacturer’s instructions. Finally, deep sequencing was performed using the Illumina Hiseq 2000 (Illumina) according to the protocol of the manufacturer.
Sequencing data analysis
The raw data were retrieved by Illumina’s Sequencing Control Studio software version 2.8 (SCS v 2.8, LC-Bio, Houston, TX, USA), and extracted from the image files generated by Illumina Genome Analyzer Pipeline software (LC Sciences, Hangzhou, China). The clean reads (18–25 nt) were collected after removing junk sequences, error length sequences (< 18 nt and > 25 nt), and sequences mapped to RFam (http://rfam.janelia.org), RepBase (http://www.girinst.org/repbase) and mRNA (http://www.ncbi.nlm.nih.gov/). Additionally, the unique clean reads mapped to miRNA precursors found using a BLAST algorithm-based search against the miRBase v21.0 (htt://mirbase.org) were identified as known miRNAs, allowing up to two mismatches during the alignment. Others were mapped to genomes of B. macrocarpa [28] or reference genomes of B. vulgaris in NCBI to identify their precursors. Those with precursors that met the criteria (MFEI > 0.9) were marked as potential miRNAs. RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi) was used to investigate their secondary hairpin structures. Other criteria are listed in S1 Table.
Phylogenetic analysis
The MIR160 (precursor of miR160) phylogenetic tree with 34 MIR160s from 29 plant species was constructed by MEGA5.0 using the Maximum Likelihood method according to the manual [50]. These validated MIR160s from various plant species were adapted as reported [9, 10, 33], and their sequences were downloaded from miRBase v21.0. The phylogenetic distribution of the 29 plant species was analyzed by the common tree tool in NCBI (https://www.ncbi.nlm.nih.gov/Taxonomy/CommonTree/). The three-letter codes of the full species names are listed.
Microarray analysis
The miRNA microarray experiment was performed according to the protocol provided by LC Sciences (LC Sciences, Hangzhou, China). The small RNAs (< 300-nt) were size-fractionated using YM-100 microcon centrifugal filters (Millipore) from 2 μg of total RNA samples. These small RNAs were 3′ extended with a poly(A) tail using poly(A) polymerase and then an oligonucleotide tag was ligated to the poly(A) tail for later fluorescent dye staining. Antisense detection probes for each of the 158 chosen identified miRNAs were made by in situ synthesis using photo-generated reagent chemistry. Hybridization was performed overnight on a μParaflo microfluidic chip using a micro-circulation pump (Atactic Technologies, Houston, TX, USA).
Each probe was used three times, arranged in different places on one chip, and three biological replicates were performed. Tag-specific Cy3 and Cy5 dyes were used to label virus-free and BNYVV-infected RNA samples, respectively, and were exchanged during the double-color microarray assays. After RNA hybridization, tag-conjugating Cy3 and Cy5 dyes were circulated through the microfluidic chip for dye staining. Fluorescence images were collected using a laser scanner (GenePix 4000B, Molecular Device, Union City, CA, USA) and digitized using Array-Pro image analysis software (Media Cybernetics, Silver Spring, MD, USA). Background signals were subtracted from the original data followed by normalization of the signals using the LOWESS filter (Locally-weighted Regression).
Degradome sequencing and analysis
One degradome library was constructed from leaves of BNYVV-infected (VL) plants based on published methods [34]. The cDNA library was sequenced on an Illumina HiSeq 2000 (LC Sciences, Hangzhou, China). CleaveLand 3.0 pipeline and ACGT301-DGEv1.0 program (LC Sciences) were employed for data analysis. The alignment scores ≤ 4 were used as the criteria. Furthermore, based on the abundance of degradome tags at the target sites, the miRNA targets were classified into 5 categories (0, 1, 2, 3, and 4) according to previously described methods [35].
Reverse-transcription quantitative real-time PCR (RT-qPCR)
To validate the small RNA sequencing results, the relative expression levels of 11 selected miRNAs from BNYVV-infected or mock-inoculated leaves were validated by stem-loop RT-qPCR. Primers for the stem–loop RT-qPCR were designed using methods as previously described [51]. Total RNA (3 μg) was used in the reverse transcription reaction (30 μl). qPCR was performed in 96-well plates using the CFX96 real-time PCR detection system (Bio-Rad, Hercules, CA, USA) with the following temperature program: 95°C for 15 s, followed by 40 cycles of 95°C for 15 s, and then annealing at 60°C for 30 s. Each reaction mixture consisted of 1 μl cDNA, 7 μl SoFast EyaGreen Supermix (Bio-Rad, Hercules, CA, USA), 1.5 μl (3 pmol/μl) of both forward and reverse primers, and 3 μl PCR-grade water [28].
The transcript levels of several target genes of six selected miRNAs were also assayed by RT-qPCR. Reverse transcription of RNA was conducted by oligo (dT20). qPCR was conducted as described above. PP2A was used as an internal control. All reactions were performed using three biological replicates for each treatment, and each biological sample of three pooled plants was evaluated with three technical replicates. All primers used in this study are listed in S6 Table.
Supporting information
S1 Table. The expression levels of miRNA in B. macrocarpa (reads > 3).
https://doi.org/10.1371/journal.pone.0186500.s001
(XLSX)
S2 Table. The list of known miRNAs-3p and -5p.
https://doi.org/10.1371/journal.pone.0186500.s002
(XLSX)
S3 Table. The list of potential novel miRNAs-3p and -5p.
https://doi.org/10.1371/journal.pone.0186500.s003
(XLSX)
S4 Table. Differential expression profiles of miRNAs between mock-, and BNYVV-infected plants (p < 0.05, Fold change > 2, reads more than 20).
https://doi.org/10.1371/journal.pone.0186500.s004
(XLSX)
S5 Table. The results of microarray analysis of 158 miRNAs compared with small RNA deep sequencing.
https://doi.org/10.1371/journal.pone.0186500.s005
(XLSX)
S1 File. The secondary structure of bma-miR160a
https://doi.org/10.1371/journal.pone.0186500.s007
(PDF)
S2 File. The secondary structure of miRn14-3p and miRn14-5p.
https://doi.org/10.1371/journal.pone.0186500.s008
(PDF)
S3 File. The secondary structure of miRn23-3p and miR23-5p.
https://doi.org/10.1371/journal.pone.0186500.s009
(PDF)
S4 File. The secondary structure of miRn31-3p and miRn31-5p.
https://doi.org/10.1371/journal.pone.0186500.s010
(PDF)
S5 File. The secondary structure of miRn33-3p and miRn33-5p.
https://doi.org/10.1371/journal.pone.0186500.s011
(PDF)
S6 File. Heatmap of differentially expressed miRNAs (P-value < 0.05).
https://doi.org/10.1371/journal.pone.0186500.s012
(PDF)
S1 Fig. The western blot analysis of BNYVV CP expression in systemic leaves.
Mock indicates healthy B. macrocarpa inoculated with buffer only.
https://doi.org/10.1371/journal.pone.0186500.s015
(TIFF)
Acknowledgments
We thank Drs. Xian-Bin Wang, Zong-Ying Zhang, and Qian Wang for their useful comments on this study.
References
- 1. Griffiths-Jones S, Saini HK, van Dongen S, Enright AJ. miRBase: tools for microRNA genomics. Nucleic Acids Research. 2008; 36: 154–158.
- 2. Reinhart BJ, Weinstein EG, Rhoades MW, Bartel B, Bartel DP. MicroRNAs in plants. Genes & Development. 2002; 16(13): 1616–1626.
- 3. Xie ZX, Allen E, Fahlgren N, Calamar A, Givan SA, Carrington JC. Expression of Arabidopsis MIRNA genes. Plant Physiology. 2005; 138(4): 2145–2154. pmid:16040653
- 4. Kim YJ, Zheng BL, Yu Y, Won SY, Mo BX, Chen XM. The role of Mediator in small and long noncoding RNA production in Arabidopsis thaliana. EMBO Journal. 2011; 30(5): 814–822. pmid:21252857
- 5. Zhang BH, Wang QL. MicroRNA-based biotechnology for plant improvement. J Cell Physiol. 2015; 230(1): 1–15. pmid:24909308
- 6. Rogers K, Chen XM. Biogenesis, turnover, and mode of action of plant microRNAs. Plant Cell. 2013; 25(7): 2383–2399. pmid:23881412
- 7. Devers EA, Branscheid A, May P, Krajinski F. Stars and symbiosis: microRNA- and microRNA*-mediated transcript cleavage involved in arbuscular mycorrhizal symbiosis. Plant Physiol. 2011; 156(4): 1990–2010. pmid:21571671
- 8. Zhang XM, Zhao HW, Gao S, Wang WC, Katiyar-Agarwal S, Huang HD, et al. Arabidopsis Argonaute 2 regulates innate immunity via miRNA393*-mediated silencing of a Golgi-localized SNARE gene, MEMB12. Molecular Cell. 2011; 42(3): 356–366. pmid:21549312
- 9. Taylor RS, Tarver JE, Hiscock SJ, Donoghue PCJ. Evolutionary history of plant microRNAs. Trends in Plant Science. 2014; 19(3): 175–182. pmid:24405820
- 10. Chavez Montes RA, de Fatima Rosas-Cardenas F, De Paoli E, Accerbi M, Rymarquis LA, Mahalingam G, et al. Sample sequencing of vascular plants demonstrates widespread conservation and divergence of microRNAs. Nat Commun. 2014; 5: 3722. pmid:24759728
- 11. Cuperus JT, Fahlgren N, Carrington JC. Evolution and functional diversification of MIRNA genes. Plant Cell. 2011; 23(2): 431–442. pmid:21317375
- 12. Gu WJ, Wang XF, Zhai CY, Xie XY, Zhou T, et al. Selection on synonymous sites for increased accessibility around miRNA binding sites in plants. Molecular Biology and Evolution. 2012; 29(10): 3037–3044. pmid:22490819
- 13. Voinnet O. Origin, Biogenesis, and activity of plant microRNAs. Cell. 2009; 136(4): 669–687. pmid:19239888
- 14. Yang J, Zhang F, Li J, Chen JP, Zhang HM. Integrative analysis of the microRNAome and transcriptome illuminates the response of susceptible rice plants to Rice Stripe Virus. PLoS One. 2016; 11 (1): e0146946. pmid:26799317
- 15. Zhang C, Ding ZM, Wu KC, Yang L, Li Y, Yang Z, et al. Suppression of Jasmonic acid-mediated defense by viral-inducible microRNA319 facilitates virus infection in rice. Molecular Plant. 2016; 9(9): 1302–1314. pmid:27381440
- 16. Du Z, Chen A, Chen W, Westwood JH, Baulcombe DC, Carr JP. Using a viral vector to reveal the role of microRNA159 in disease symptom induction by a severe strain of Cucumber mosaic virus. Plant Physiol. 2014; 164(3): 1378–88. pmid:24492335
- 17. Jay F, Wang Y, Yu A, Taconnat L, Pelletier S, Colot V, et al. Misregulation of AUXIN RESPONSE FACTOR 8 underlies the developmental abnormalities caused by three distinct viral silencing suppressors in Arabidopsis. PLoS Pathogens. 2011; 7(5).
- 18. Tamada T, Baba T. Beet necrotic yellow vein virus from rizomania-affected sugar beet in Japan. Ann Phytopathol Soc Jpn. 1973; 39: 325–332.
- 19. Tamada T, Shirako Y, Abe H, Saito M, Kiguchi T, Harada T. Production and pathogenicity of isolates of Beet necrotic yellow vein virus with different numbers of RNA components. J Gen Virol. 1989; 70: 3399–3409.
- 20. Richards KE, Tamada T. Mapping functions on the multipartite genome of Beet Necrotic Yellow Vein Virus. Annual Review of Phytopathology. 1992; 30: 291–313.
- 21. Tamada T, Uchino H, Kusume T, Saito M. RNA 3 deletion mutants of beet necrotic yellow vein virus do not cause rhizomania disease in sugar beets. Phytopathology. 1999; 89(11): 1000–1006. pmid:18944654
- 22. Rahim MD, Andika IB, Han C, Kondo H, Tamada T. RNA4-encoded p31 of Beet necrotic yellow vein virus is involved in efficient vector transmission, symptom severity and silencing suppression in roots. Journal of General Virology. 2007; 88: 1611–1619. pmid:17412994
- 23. McGrann GRD, Grimmer MK, Mutasa-Goettgens ES, Stevens M. Progress towards the understanding and control of sugar beet rhizomania disease. Molecular Plant Pathology, 2009; 10(1): 129–141. pmid:19161359
- 24. Chiba S, Kondo H, Miyanishi M, Andika IB, Han CG, Tamada T. The evolutionary history of Beet necrotic yellow vein virus deduced from genetic variation, geographical origin and spread, and the breaking of host resistance. Molecular Plant-Microbe Interactions. 2011; 24(2): 207–218. pmid:20977309
- 25. Thiel H, Hleibieh K, Gilmer D, Varrelmann M. The P25 pathogenicity factor of Beet necrotic yellow vein virus targets the sugar beet 26S proteasome involved in the induction of a hypersensitive resistance response via interaction with an F-box protein. Molecular Plant-Microbe Interactions. 2012; 25(8): 1058–1072. pmid:22512382
- 26. Thiel H, Varrelmann M. Identification of Beet necrotic yellow vein virus P25 pathogenicity factor-interacting sugar beet proteins that represent putative virus targets or components of plant resistance. Molecular Plant-Microbe Interactions. 2009; 22(8): 999–1010. pmid:19589075
- 27. Peltier C, Schmidlin L, Klein E, Taconnat L, Prinsen E, Erhardt M, et al. Expression of the Beet necrotic yellow vein virus p25 protein induces hormonal changes and a root branching phenotype in Arabidopsis thaliana. Transgenic Res. 2011; 20(3): 443–66. pmid:20602166
- 28. Fan HY, Zhang YL, Sun HW, Liu JY, Wang Y, Wang XB, et al. Transcriptome analysis of Beta macrocarpa and identification of differentially expressed transcripts in response to Beet Necrotic Yellow Vein Virus infection. PLoS One. 2015; 10(7): e0132277. pmid:26196682
- 29. Baksa I. Nagy T. Barta E. Havelda Z. Varallyay E. Silhavy D, et al. Identification of Nicotiana benthamiana microRNAs and their targets using high throughput sequencing and degradome analysis. BMC Genomics. 2015; 16.
- 30. Wang TZ, Chen L, Zhao MG, Tian QY, Zhang WH. Identification of drought-responsive microRNAs in Medicago truncatula by genome-wide high-throughput sequencing. BMC Genomics. 2011; 12.
- 31. Saminathan T, Bodunrin A, Singh NV, Devarajan R, Nimmakayala P, Jeff M, et al. Genome-wide identification of microRNAs in pomegranate (Punica granatum L.) by high-throughput sequencing. BMC Plant Biol. 2016; 16(1): 122. pmid:27230657
- 32. Fahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, Cumbie JS, et al. High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes. PLoS One. 2007; 2(2).
- 33. Yin HF. Fan ZQ. Li XL. Wang JY. Liu WX. Wu B, et al. Phylogenetic tree-informed microRNAome analysis uncovers conserved and lineage-specific miRNAs in Camellia during floral organ development. J Exp Bot. 2016; 67(9): 2641–53. pmid:26951373
- 34. Addo-Quaye C, Eshoo TW, Bartel DP, Axtell MJ. Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Current Biology. 2008; 18: 758–762. pmid:18472421
- 35. Addo-Quaye C, Miller W, Axtell MJ. CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics. 2009; 25: 130–131. pmid:19017659
- 36. Varallyay E, Valoczi A, Agyi A, Burgyan J, Havelda Z. Plant virus-mediated induction of miR168 is associated with repression of ARGONAUTE1 accumulation. EMBO J. 2010; 29(20): 3507–3519. pmid:20823831
- 37. Gang L, Yang F, Yu D. Microrna395 mediates regulation of sulfate accumulation and allocation in Arabidopsis thaliana. Plant Journal. 2010; 62(6): 1046–1057. pmid:20374528
- 38. Yu S, Galvao VC, Zhang YC, Horrer D, Zhang TQ, Hao YH, et al. Gibberellin regulates the Arabidopsis floral transition through miR156-targeted SQUAMOSA PROMOTER BINDING-LIKE transcription factors. Plant Cell. 2012; 24(8): 3320–3332. pmid:22942378
- 39. Eviatar-Ribak T, Shalit-Kaneh A, Chappell-Maor L, Amsellem Z, Eshed Y, Lifschitz E. A cytokinin-activating enzyme promotes tuber formation in Tomato. Current Biology. 2013; 23(12): 1057–1064. pmid:23746638
- 40. Xie K, Wu C, Xiong L. Genomic organization, differential expression, and interaction of SQUAMOSA promoter-binding-like transcription factors and microRNA156 in rice. Plant Physiol. 2006; 142(1): 280–293. pmid:16861571
- 41. Fan G, Cao X, Niu S, Deng M, Zhao Z, Dong Y. Transcriptome, microRNA, and degradome analyses of the gene expression of Paulownia with phytoplamsa. BMC Genomics. 2015; 16: 896. pmid:26537848
- 42. Gai YP, Li YQ, Guo FY, Yuan CZ, Mo YY, Zhang HL, et al. Analysis of phytoplasma-responsive sRNAs provide insight into the pathogenic mechanisms of mulberry yellow dwarf disease. Scientific Reports. 2014; 4.
- 43. Poethig RS. Small RNAs and developmental timing in plants. Curr Opin Genet Dev. 2009; 19(4): 374–378. pmid:19703647
- 44. Fornara F and Coupland G. Plant phase transitions make a SPLash. Cell. 2009; 138(4): 625–627. pmid:19703391
- 45. Xu MY, Zhang L, Li WW, Hu XL, Wang MB, Fan YL, et al. Stress-induced early flowering is mediated by miR169 in Arabidopsis thaliana. J Exp Bot, 2014; 65(1): p. 89–101. pmid:24336445
- 46. Li SB, Xie ZZ, Hu CG, Zhang JZ. A review of Auxin Response Factors (ARFs) in plants. Front Plant Sci. 2016; 7: 47. pmid:26870066
- 47. Si-Ammour A, Windels D, Arn-Bouldoires E, Kutter C, Ailhas J, Meins F, et al. miR393 and secondary siRNAs regulate expression of the TIR1/AFB2 auxin receptor clade and auxin-related development of Arabidopsis leaves. Plant Physiology. 2011; 157(2): 683–691. pmid:21828251
- 48. Schommer C, Palatnik JF, Aggarwal P, Chetelat A, Cubas P, Farmer EE, et al. Control of jasmonate biosynthesis and senescence by miR319 targets. PLoS Biology. 2008; 6(9): 1991–2001.
- 49. Wang Y, Fan HY, Wang XB, Li M., Han CG. Li DW, et al. Detection and characterization of spontaneous internal deletion mutants of Beet Necrotic yellow vein virus RNA3 from systemic host Nicotiana benthamiana. Virology Journal. 2011; 8.
- 50. Tamura K, Peterson D, Peterson N, Stecher G, Ne M, Kumar S. Mega5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Molecular Biology & Evolution, 2011; 28(10): 2731–2739.
- 51. Chen C, Ridzon DA, Broomer AJ, Zhou Z, Lee DH, Nguyen JT, et al. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucleic acids research. 2005; 33: e179. pmid:16314309