Identification of microRNAs regulated by tobacco curly shoot virus co-infection with its betasatellite in Nicotiana benthamiana

Background MicroRNAs (miRNAs) are a class of 21–24 nucleotide endogenous non-coding small RNAs that play important roles in plant development and defense responses to biotic and abiotic stresses. Tobacco curly shoot virus (TbCSV) is a monopartite begomovirus, cause leaf curling and plant stunting symptoms in many Solanaceae plants. The betasatellite of TbCSV (TbCSB) induces more severe symptoms and enhances virus accumulation when co-infect the plants with TbCSV. Methods In this study, miRNAs regulated by TbCSV and TbCSB co-infection in Nicotiana benthamiana were characterized using high-throughput sequencing technology. Results Small RNA sequencing analysis revealed that a total of 13 known miRNAs and 42 novel miRNAs were differentially expressed in TbCSV and TbCSB co-infected N. benthamiana plants. Several potential miRNA-targeted genes were identified through data mining and were involved in both catalytic and metabolic processes, in addition to plant defense mechanisms against virus infections according to Gene Ontology (GO) analyses. In addition, the expressions of several differentially expressed miRNAs and their miRNA-targeted gene were validated through quantitative real time polymerase chain reaction (qRT-PCR) approach. Conclusions A large number of miRNAs are identified, and their target genes, functional annotations also have been explored. Our results provide the information on N. benthamiana miRNAs and would be useful to further understand miRNA regulatory mechanisms after TbCSV and TbCSB co-infection.


Background
MicroRNAs (miRNAs) are a large group of small, noncoding RNAs that are 21-24 nucleotides (nt) long and regulate gene expression at the post-transcriptional level [1,2]. In plants, miRNA is initially produced by RNA polymerase II (Pol II) or III (Pol III) and is known as the primary miRNA (pri-miRNA). Then the pri-miRNA is processed by a Dicer-like (DCL) enzyme, which produces an miRNA precursor (pre-miRNA), which contains both a miRNA:miRNA* intermediate duplex and a stem-looplike structure [3]. After excision, the mature miRNA is bound to an argonaute (AGO) protein to form a RNA-induced silencing complex (RISC) [4]. Finally, the RISC complex recognizes and binds with the target RNA transcript, in a sequence-specific manner, and then cleaves the transcript or represses the translation process [2].
In addition, many studies have indicated that miRNAs are highly conserved in eukaryotes, and virus infections in plants often alter miRNA expressions. For example, tobacco mosaic virus (TMV) infection in Nicotiana tabacum changes the expressions of several specific miR-NAs [17]. Infection of cotton plants with cotton leaf roll dwarf virus (CLRDV) affects miRNA expression, downregulating the expression of specific miRNAs and causing disease symptoms in cotton leaves [18]. In A. thaliana and N. benthamiana plants, the expressions of miR168 and the Argonaute 1 gene are upregulated after infection with several plant viruses [19]. The expression of miR159 is upregulated by cucumber green mottle mosaic virus (CGMMV) infection cucumber at 10 post-inoculation [20]. In rice, downregulation of osa-miR171b expression in plants infected with rice stripe virus (RSV) enhances RSV disease symptoms [21]. The expression of miR319 in rice is increased after infection with rice ragged stunt virus (RRSV), and upregulation of miR319 promotes RRSV infection and disease symptoms by inhibiting JA-mediated host defense [22]. Amin et al. (2011) reported that begomovirus infection upregulates the accumulation of miR-NAs controlling plant development [23].
In this study, a miRNA-sequencing approach was used to identify N. benthamiana miRNAs regulated by TbCSV/ TbCSB co-infection. Then we predicted the target genes of the differentially expressed miRNAs. To investigate the interactions between the differentially expressed miRNAs and their target genes during TbCSV/TbCSB co-infection in N. benthamiana, the expression levels of several identified miRNAs together with their target genes were analyzed using qRT-PCR. The results not only shed light on the possible roles of miRNAs in N. benthamiana development and physiology, but also their possible roles in N. benthamiana resistance to TbCSV/TbCSB co-infection.

Plant growth and virus inoculation
N. benthamiana plants were grown inside a greenhouse at 24°C and a 16/8 h (light/dark) photoperiod. Infectious clones of TbCSV isolate Y35 and its betasatellite (TbCSB) were individually introduced to Agrobacterium tumefaciens strain EHA105 [33]. After overnight culturing and centrifugation, agrobacterium pellets were resuspended in infiltration buffer (10 mM MES, 10 mM MgCl 2 and 200 μM acetosyringone in sterile water) until they reached an OD 600 value of 2, followed by 3 h incubation at room temperature. The two agrobacterium cultures were mixed 1:1 (v/v) and then co-infiltrated into the leaves of N. benthamiana plants at the six to eight leaf stage as described previously [35]. Plants infiltrated with the infiltration buffer only were used as noninfected controls.

Sample preparation and total RNA isolation
At 20 days-post co-infiltration (dpi), systemic leaves of three TbCSV/TbCSB-infected and three non-infected control plants were harvested. Total RNA was isolated from leaf samples using the TRIzol reagent, following the manufacturer's instructions (Invitrogen, Waltham, MA, USA).

Small RNA library construction and Illumina sequencing
The integrity and concentration of RNA samples were checked using a NanoDrop spectrophotometer and an Agilent 2100 Bioanalyzer, following the manufacturer's instructions (Agilent, Santa Clara, CA, USA). Highquality samples were identified and used to construct small RNA libraries followed by Illumina sequencing using an IlluminaHiseq™2500 instrument. Small RNA libraries and Illumina sequencing were performed by the Novogene Bioinformatics Technology Company, Beijing, China.

Analyses of Illumina sequencing date
Raw reads were processed to remove both adaptors and low-quality reads. Both the commercial Gene-Chip® Tomato Genome Array and the miRBase (http://www.mirbase.org/) were searched to identify known miRNAs in the two small RNA libraries (i.e., TbCSV/TbCSB and mock). The reads matched to the N. benthamiana genome shotgun-sequence assemblies were kept for further identifications.

Prediction of miRNA-targeted genes and gene function analyses
Potential miRNA-targeted genes were predicted using the psRNATarget program (http://plantgrn.noble.org/ psRNATarget). The rules used for the predictions were as described by Schwab et al. [36]. To explore the possible functions of the predicted target genes, Gene Ontology (GO) analyses were performed as described previously [37]. The results were split into three different categories: Biological process, Molecular functions and Cellular components. Through this process, all potential target genes were mapped to GO terms described in the database (http://www.geneontology.org).
Validations of miRNA and target gene expression qRT-PCR was performed to validate the expressions of miRNAs together with their target genes. Total RNA was isolated at 20 dpi from N. benthamiana leaves either infected with TbCSV/TbCSB or not, using TRIzol reagent. For miRNAs, the cDNA was synthesized by RT using the PrimeScript RT reagent Kit with gDNA Eraser with a special stem-loop RT primer according to the manufacturer's instructions (TaKaRa, Dalian, China). The cDNA products were used as templates for realtime PCR analyses. The primers used for qRT-PCR assays are listed in (Additional file 1 Table S1). In addition, we also used qRT-PCR to analyze the expression of target genes, and the primers are also listed in Table S1. Purified total RNA was used as a template and reversetranscribed using the PrimeScript RT reagent Kit (TaKaRa, Dalian, China) to obtain cDNA. The reactions consisted of incubation in 96-well plates at 98°C for 2 min, followed by 40 cycles of 98°C for 10 s, 60°C for 10 s, and 68°C for 30 s. The expression level of the N. benthamiana Ubiquitin C gene (UBC) was used as the reference [38,39]. Three biological samples were used for each treatment and each biological sample had a further three technical replicates during qPCR. The qPCR results were calculated using the 2 -ΔΔCT method [38].

Small RNA library construction
Small RNA libraries were constructed using total RNA from N. benthamiana leaves either infected with TbCSV/TbCSB or not. A total of 16,753,586 and 17, 822,708 raw reads were generated via Illumina sequencing to represent these samples. After removing raw reads without 3′ adapters, with 5′ added adapters, with more than 10% unidentified nucleotides (nt), with poly (A/T/G/C) stretches, and that were less than 18 nt or more than 30 nt, a total of 13,896,237 and 15,227,477 clean reads were obtained from infected and noninfected libraries, respectively (Table 1) The resulting clean small RNA reads belonged to different categories, including exon sense and antisense, intron sense and antisense, rRNA, tRNA, snRNA, snoRNA, RNA repeats, natural antisense transcripts (NAT), trans-acting siRNA (TAS), and other unannotated reads. Among these, there were 404,032 (4.01%) miRNA tags from the infected library and 408,746 (3.22%) from the non-infected library ( Table 2).
When all of the reads (18 to 30 nt) were analyzed, those with 24 nt were the most abundant. Of these, 6, 119,284 (44.04%) were found in the infected library and 7,109,221 (46.69%) were found in the non-infected library. The second most abundant reads (1,760,916 reads or 11.56%) in the non-infected library had 21 nt, and that (2,079,345 reads or 14.96%) in the infected library had 22 nt (Fig. 2 a). When only unique reads were considered, 24 nt reads were the most abundant class with a total of 3,538,301 reads (62.34%) in the infected library and 3,896,777 reads (65.14%) in the noninfected library. The second most abundant unique read (510,867 reads or 9.0%) in the infected library had  22 nt while that (438,923 reads or 7.34%) in the noninfected library had 23 nt (Fig. 2b).

Identification of known miRNAs
To identify known miRNAs in the two libraries, all small RNA reads were used to blast search the miRBase site with known mature plant miRNAs and then the N. benthamiana genome database. Through this approach, a total of 349,937 small RNA reads from the infected library and 346,265 reads from the non-infected library were mapped to the Solanum lycopersicum genome.
Excluding miRNAs expressed at extremely low levels, 40 miRNAs in 21 known miRNA families were identified. The number of small RNA reads matched to known miRNA families are summarized in Figs. 3 and 4. The expression levels of miRNAs changed slightly after infection, as shown in Table 3. There were 13 differentially expressed miRNAs within the two libraries, with 2 downregulated miRNAs and 11 upregulated miRNAs after TbCSV/TbCSB co-infection.

Identification of novel candidate miRNAs
To predict hairpin-like structures in the identified miRNA precursors and identify the corresponding miR-NAs that could be used to further identify novel miR-NAs, we utilized the miREvo and miRDeep2 software as described [40,41]. We also used the mfold web server (http://unafold.rna.albany.edu) to predict the secondary structures and the minimum free energy of the annotated small RNA tags that were mapped to the N. benthamiana genome, as described previously [42]. After removing miRNAs with extremely low expressions, a total of 42 miRNAs were found to show differential   Table 3, Additional file 2 Figure  S1). Among the predicted novel miRNA candidates, the base bias in the first position showed that the majority of these novel miRNA candidates started with a 5′ uridine (U) as shown in Fig. S2a and Fig. S2b. Furthermore, novel miRNA candidate nucleotide bias at each position were also analyzed (Additional file 3 Figure S2c and Figure S2d). In addition, we also found among the 42 differentially-expressed novel miRNAs, twenty-five had complementary miRNA*s, with precursor lengths ranging from 43 to 295 nt and predicted minimal folding energy (MFE) ranging from − 14.2 to − 130.7 kcal/mol.    Table  S2). All of the novel miRNAs' loci, pre-miRNA sequences and structures, and reads from deep sequencing were shown in Additional file 4. This is in agreement with published criteria for novel miRNA [43,44], and suggests that these candidate miRNAs are most likely to be new miRNA family members in N. benthamiana.

Prediction of potential miRNA-targeted genes
Numerous genes are responsive to virus infections and to differentially expressed miRNAs. Our results indicate that many potential miRNA-targeted genes encode transcription or non-transcription factor proteins, which are important for physiological processes. To explore the regulatory functions of the identified miRNAs in the infected library, potential target genes of nine conserved and four novel miRNAs were predicted by GO analyses. The GO annotated terms Biological process, Cellular components, and Molecular function were further analyzed to determine genes that could potentially be targeted by the identified miRNAs (Fig. 5). For the

Validations of miRNA and target gene expressions by qRT-PCR
To validate the high-throughput sequencing results, nine known and four novel miRNAs that showed differential expression between the two libraries were selected and analyzed for expression by stem-loop qRT-PCR (Table 4 and Fig. 6). The PCR primers are listed in Table S1. The expressions of miR156d-5p, miR169c, miR4376, miR156a, miR1919c-5p, miR159, novel-121, novel-71, and miR482a in the infected library were all upregulated, whereas the expressions of miR171b, miR164a-5p, novel-94, and novel-70 were all downregulated. To examine if the expressions of TbCSV/TbCSB infection-regulated miRNAs could influence the expressions of their target genes, we analyzed the predicted target genes through qRT-PCR. The expressions of squamosa promoter-binding-like protein (targeted by miR156a), disease resistance protein (targeted by miR482a), nuclear transcription factor Y subunit (targeted by miR169c), calcium-transporting ATPase (targeted by miR4376), conserved hypothetical protein (targeted by miR4376), MYB-like transcription factor (targeted by miR159), heavy metal transport (targeted by novel-121), and N-acetyltransferase (targeted by novel-71) were all downregulated after TbCSV/TbCSB infection. By contrast, GRAS family transcription factor (targeted by miR171), NAC domain-containing protein (targeted by miR164a), transcription factor (targeted by novel-70), and 1-aminocyclopropane-1-carboxylate oxidase (targeted by novel-94) were all upregulated after TbCSV/TbCSB coinfection (Table 4 and Fig. 7).

Discussion
High-throughput sequencing technology has been used extensively in small RNA research [45]. More and more studies have illustrated that various virus infections in plants often alter miRNA expressions [17][18][19][20]46]. A large number of miRNAs have been identified in plants and the functions of many miRNAs have also been investigated [47,48]. To better understand the roles of the N. benthamiana miRNAs in host resistance to TbCSV/ TbCSB co-infection, in this study, two libraries were constructed, using total RNA from N. benthamiana plants either infected with TbCSV/TbCSB or not.
The most abundant small RNA reads in the two libraries were those with 24 nt. In addition, the 24 nt small RNA reads were the predominant unique small RNA reads. This finding agrees with earlier reports which have shown that 24 nt small RNAs are more abundant in several other diseased plants, such as tomato plants infected with Phytophthora infestans, tomato infected with cucumber mosaic virus (CMV) and wheat plants infected with powdery mildew pathogen [13,49,50]. The length distribution of small RNA reads may reflect their compositions [51]. We found more 24 nt small RNAs in the non-infected library than in the infected library. By contrast, 21 and 22 nt miRNAs were more abundant in the infected library. Our results also indicate that the expression profiles of miRNAs were significantly altered after TbCSV/TbCSB co-infection in N. benthamiana, and the differentially regulated expressions of miRNAs Fig. 7 Relative expressions of miRNA-targeted genes. Relative expression of each gene was determined through qRT-PCR. a Relative expressions of the known miRNA-targeted genes. b Relative expressions of the novel miRNA-targeted genes. The x-axis shows the names of the miRNAtargeted genes analyzed in this study. The y-axis shows the Log 2 ratio between the expression values from a TbCSV/TbCSB-infected sample versus its Mock sample. Three biological replicates were analyzed for each miRNA-targeted gene through qRT-PCR. Expression level of N. benthamiana UBC gene was used as the reference gene during qRT-PCR assays. Asterisks indicate statistically significant differences compared with Mock, "*" indicate significant difference (0.01 ≤ p < 0.05), "**" indicate extremely significant difference (P < 0.01) suggested that miRNAs play important roles during TbCSV/TbCSB co-infection.
In all, 13 known and 42 potentially novel miRNAs were differentially regulated by TbCSV/TbCSB coinfection. To better understand the relative abundances of miRNAs in the two libraries, we analyzed sequence frequencies and used them as indexes. When the two libraries were compared, the normalized reads varied from about 2 (novel-96) to 63,693 (miR482a) in the infected library and from 0 (novel-131) to 28,999 (miR482a) in the non-infected control library, indicating significant variation in the relative abundances of different miRNA sequences. This finding was later confirmed through qRT-PCR analyses using several selected differentially expressed miRNAs. miRNAs may regulate host defenses against pathogens, including viruses, by suppressing pathogen multiplication at the post-transcriptional level [13,52]. Several stress-responsive miRNAs (e.g., miR168, miR169, and miR482) have been reported to target transcription factors controlling host resistance to virus infection [53][54][55]. In N. benthamiana, virus infection may regulate the expression of miR168 to alleviate the anti-viral function of AGO1 protein [53]. In our study, N. benthamiana miR168 was found to be responsive to TbCSV/TbCSB co-infection. In addition, the expression of miR169 was upregulated after co-infection. In a previous study, rice miR169 was found to negatively regulate rice immunity against Magnaporthe oryzae infection by differentially repressing its target genes [54]. Studies also showed that rice miR164 plays an important role in rice resistance to southern rice black-streaked dwarf virus (SRBSDV) infection as well as rice resistance to drought stresses by differentially regulating its target genes [56,57]. In addition, miR482 can regulate the expression of NBS-LRR defense genes during fungal pathogen infection in cotton [55]. In our study, when the expressions of miR164a and miR482 in the two libraries were compared, the normalized miRNA164a and miR482 reads were 3624 and 28,999 in the noninfected control library, and 1633 and 63,693 in the infected library, respectively, suggesting that these two miRNAs may play roles in N. benthamiana resistance to TbCSV/TbCSB co-infection. The further studies will be continued to unravel the functions of these miRNAs.

Conclusion
In this study, miRNAs regulated by TbCSV and TbCSB co-infection in N. benthamiana were characterized using high-throughput sequencing technology, and some miR-NAs involved in plant defense system were found to be significantly regulated after TbCSV and TbCSB infection. The molecular functions of these miRNAs in N.
benthamiana resistance to TbCSV/TbCSB co-infection may require further investigation. Nonetheless, our results improve knowledge of the infection of TbCSV/ TbCSB in host plants while also providing additional information for the development of management strategies for TbCSV/TbCSB infection in the future.