Characterization of microRNAs of Beta macrocarpa and their responses to Beet necrotic yellow vein virus infection

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.


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][21][22][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 miR-NAs 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 virusfree 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.
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

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 ).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.
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 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 miR-NAs 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).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 (|log 2 (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.
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].
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][40][41][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 BNYVVinfected 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.

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 K 2 HPO 4 , 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

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

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 0 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 doublecolor 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 miR-NAs 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 (dT 20 ).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.

Fig 3 .
Fig 3.The distribution of miRNAs of B. macrocarpa within the phylogenetic perspective of plant miRNAs.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.https://doi.org/10.1371/journal.pone.0186500.g003

Fig 4 .
Fig 4. The phylogenetic distribution of the 29 plant species analyzed in this study.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.https://doi.org/10.1371/journal.pone.0186500.g004

Fig 7 .
Fig 7. Validation of the relative expression levels of selected miRNAs and their targets by RT-qPCR.(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.https://doi.org/10.1371/journal.pone.0186500.g007

Table 1 . Summary of small RNA sequencing from the BNYVV-infected and virus-free library of B. macrocarpa. category a L VL Total Unique Total Unique Raw reads
Table, S5 Table in bold font).Additionally, 10 of 15 other miRNA families were tested again by microarray analysis (S5 Table in bold font