The First Report of miRNAs from a Thysanopteran Insect, Thrips palmi Karny Using High-Throughput Sequencing

Thrips palmi Karny (Thysanoptera: Thripidae) is the sole vector of Watermelon bud necrosis tospovirus, where the crop loss has been estimated to be around USD 50 million annually. Chemical insecticides are of limited use in the management of T. palmi due to the thigmokinetic behaviour and development of high levels of resistance to insecticides. There is an urgent need to find out an effective futuristic management strategy, where the small RNAs especially microRNAs hold great promise as a key player in the growth and development. miRNAs are a class of short non-coding RNAs involved in regulation of gene expression either by mRNA cleavage or by translational repression. We identified and characterized a total of 77 miRNAs from T. palmi using high-throughput deep sequencing. Functional classifications of the targets for these miRNAs revealed that majority of them are involved in the regulation of transcription and translation, nucleotide binding and signal transduction. We have also validated few of these miRNAs employing stem-loop RT-PCR, qRT-PCR and Northern blot. The present study not only provides an in-depth understanding of the biological and physiological roles of miRNAs in governing gene expression but may also lead as an invaluable tool for the management of thysanopteran insects in the future.


Introduction
MicroRNAs (miRNAs) are a family of small (~18-25 nucleotides (nts)), endogenously initiated, non-coding RNAs (ncRNAs) that primarily regulate gene expression in animals, plants and protozoan in a sequence-specific manner. In mammals, approximately 60% of proteincoding gene activities are under the control of miRNAs and they regulate almost every cellular process investigated [1,2]. miRNAs can regulate gene expression either by translation repression or by degradation of mRNA through deadenylation [3]. The second to seventh nucleotides in the 5' end of the miRNA form the "seed" region that provides the most of the pairing specificity [4,5]. miRNA-mediated regulation plays a key role in cellular and developmental the subject species in this study; T. palmi is an invertebrate and not listed on the endangered species list. Thrips are ubiquitous in their natural ranges.

Insect Culture
T. palmi were reared on French bean pods (Phaseolus vulgaris cv. Arka Komal) in plastic containers (10X10 cm) at 30 ± 2°C, 80 ± 10% RH and 16:8 L:D photoperiod as described previously [36]. Total RNA was isolated from whole-body homogenates of sample mix, containing a total of 50 mg of different life stages such as eggs, larvae, pupae and adults of T. palmi using TRIzol reagent (Invitrogen, Carlsbad CA, USA).

Library preparation and sequence data generation
Samples were processed according to Illumina TruSeq™ Small RNA sample preparation guide. Size fractionated small RNA populations (18-28 nts) were extracted, purified and ligated to 3' and 5' adapters using T4 RNA Ligase (Life Technologies, Ambion, USA). Ligated products were reverse transcribed using SuperScript II (Life Technologies, Invitrogen, USA) followed by PCR amplification with 11 cycles and two size selection gels. High-throughput sequencing of the small RNA libraries was performed on Illumina Hiseq2000.

Bioinformatics analysis
The obtained sequence tags were subjected to a primary analysis in which low-quality tags, 3' and 5' adapter contaminants were discarded. The sequencing data were investigated against the Rfam (http://rfam.sanger.ac.uk/) and RepBase (http://www.girinst.org/repbase/) as references to annotate the ncRNAs namely, rRNAs, tRNAs, snRNAs, snoRNAs and repeat-associated small RNAs and degraded fragments of expressed genes (exons and introns) in the remaining sequences. All such mapped reads were removed from the dataset before further analysis. Remaining unique sequences were aligned with the miRNA sequences available in the miRBase (v21, http://www.miRBase.org/) to identify the conserved miRNAs. Novel miRNA candidates were identified by employing the miRDeep2 [37] and miRCat (http://srna-workbench.cmp. uea.ac.uk/tools/mircat/) software. Frankliniella occidentalis genome (http://www.ncbi.nlm.nih. gov/nuccore/644576459) was used as a reference to extract the potential secondary hairpin structures, employing Vienna RNAfold [38].

Homology analysis
Homology analysis was carried out with conserved miRNAs of T. palmi with the miRNAs of other organisms from the miRBase database (Release 21.0) [39]. BLASTn embedded in the miRBase database was used to compare the T. palmi miRNAs with other species, with an Evalue of 0.01 to find out more miRNA homologs. The naming of the miRNAs in this study has been done according to Griffith-Jones, et al., 2006. Since these miRNAs were predicted from T. palmi, the prefix for all miRNAs was fixed as 'tpa' . The rest of the naming convention criteria were in accordance with miRBase [40].

Prediction of miRNA targets
Target identification is crucial for understanding the biological functions of miRNAs. Unlike the plant counterparts, the imperfect complementarity of animal miRNAs with their target sequences on mRNA makes it more difficult to judge the accuracy of prediction [20]. Targets for identified miRNAs were predicted employing miRanda program [42], against the Expressed Sequence Tags (ESTs) and transcriptome (NCBI Accession: PRJNA203209) database of F. occidentalis. An alignment score [43] greater than or equal to 130 and miRNA: mRNA binding energy (Minimum Free Energy (MFE, ΔG)) less than -20 kcal/mol were considered as putative target genes. The targets were further annotated against NCBI-RefSeq invertebrate protein database and Gene Ontology (GO) terms were assigned (using Blast-2-GO) based on the annotation. The circos plot was generated using Circos [44] to visualize the interaction between miRNAs and their targets.
Validation and quantification of T. palmi miRNAs stem-loop RT-PCR. We were able to validate few of the conserved and novel microRNAs employing Stem-loop RT-PCR primers designed based on the previous reports [22]. Briefly, stem-loop RT primers bind to the 3' portion of miRNA molecules, initiating reverse transcription of the miRNAs. Later, the RT product was amplified using a miRNA specific forward primer and the universal reverse primer.
Reverse transcription quantitative PCR (qRT-PCR). In the present study, we selected differentially expressed and functionally significant 13 miRNAs (nine conserved and four novel) for qRT-PCR. Briefly, Mir-X-miRNA qRT-PCR SYBR Kit (Clontech Laboratories, Inc., USA) was used for the qRT-PCR reactions, which has a single-step, single-tube reaction to produce the first-strand cDNA, which was then specifically and quantitatively amplified using a miRNA-specific primer and SYBR Advantage qPCR chemistry. All the qRT-PCR assays were conducted according to the MIQE guidelines [45]. U6 snRNA was used as an internal control gene for normalization. qRT-PCR was performed on Light Cycler 480 (Roche, USA) using 1:20 diluted cDNAs and SYBR Advantage Premix (Clontech Laboratories, Mountain View, USA), according to the manufacturer's instructions. Assays were performed in triplicates for three independent biological experiments and the relative gene expression data were analyzed using 2 -ΔΔCT method [46]. The values of these three independent experiments were statistically analyzed using one-way ANOVA to calculate the statistical significance. miRNA northern blot. Small RNAs were isolated from pooled T. palmi (eggs, larvae, pupae and adults) employing mirVANA miRNA isolation kit (Life Technologies, USA). 5 μg of T. palmi small RNAs were resolved on 15% polyacrylamide gel containing 8M urea. RNA was electro blotted [TransBlot SD Semi-Dry Electrophoretic transfer cell (Bio-Rad, USA)] for 90 minutes at 20V onto HyBond-N+ membrane (GE Healthcare, USA) and immobilized by UV cross-linking (Stratagene, USA). The membrane was hybridized with 5' digoxigeninlabeled locked nucleic acid probes for miRNA detection (100ng/ml, Exiqon) at 37°C overnight. Later the membranes were washed twice in 2x SSC at 37°C for 15 minutes each. Digoxigenin signals were detected with DIG Northern starter kit (Roche) according to the manufacturer's instructions.

Overview of the small RNA Library
We obtained a dataset of about 14 million reads from the pooled T. palmi small RNA library (egg, larva, pupa and adult) sequenced on Illumina Next Generation Sequencing platform.
After various mapping (Table 1), the trimmed high-quality small RNA reads were employed to identify both known and novel miRNAs. Size distribution of the high-quality reads in the library (Fig 1) revealed that the peak was at 25 nts, which was also observed in the pooled library of Plutella xylostella (L.) [47]. A small portion (<5%) of our library consisted of read length of around 26 to 28 nts, which could be putative piwi-interacting RNAs (piRNAs) from T. palmi (S1 Table).

Identification of known miRNA
MiRNAs are known to be conserved among different species within a kingdom. Here, in our study, the mappable sequences were aligned to miRNA sequences from miRBase v.21.0. The analysis resulted in a total of 67 conserved miRNAs representing 54 different miRNA families (Table 2), among which the average similarity between the homologs reached 85% and few of them had a similarity to the extent of 95-100% with 1-2 nts or no difference. Analysis of the 54 miRNA families revealed that 15 were found to be exclusively present in arthropod species (Table 3), while 25 miRNA families were vertebrate specific. Seven miRNA families (miR-10, miR-100, miR-71, miR-9, miR-92, miR-15 and miR-281) were found to be highly conserved in the Animal Kingdom (Table 3) during evolution implicating their importance in regulating the gene transcripts involved in the physiological process. Among the known miRNAs, miR-281 and miR-750 were highly expressed with an expression value of 9560 and 5849 respectively ( Table 2).

Identification of miRNA-star strands
In most of the cases, once the mature miRNA strand is loaded into RISC, its star strand will be degraded soon after being exported to the cytosol. However, our analysis revealed that two T. palmi miRNA star (miRNA Ã ) strands namely, miR-6489 Ã and miR-6493 Ã were obtained from our library for their corresponding mature miRNAs (Fig 2). The expression values (Number of reads) of miR-6489 Ã were lower than that of their corresponding miRNAs, whereas, for miR-6493 Ã it was three times higher ( Table 2).

Identification of novel miRNAs
We utilized the genomic sequence assembly of F. occidentalis (Thripidae: Thysanoptera) to identify the novel miRNAs, as no genomic information is available for T. palmi. Miranalyzer pipeline identified a total of 10 novel miRNAs from T. palmi for the first time (Table 4), with their predicted precursor secondary structures (Fig 3). The complete details of the mature miR-NAs and their corresponding pre-miRNAs have been given in Table 4. The length of the novel miRNAs ranged from 21-24 nucleotides with a preference of Uracil (60%) followed by Cytosine (20%) at the 5' end. Among these ten miRNAs, five were located in the 5' arm while the other five arose from 3' arm (Table 4, Fig 3).   Abundance of novel miRNAs The novel miRNAs identified from T. palmi varied in their expression values in the library. Among the novel miRNAs, tpa-miR-N3 (4007 copies), tpa-miR-N4 (1026 copies) and tpa-miR-N7 (380 copies) had the highest abundance compared to the remaining novel miRNAs (Table 4). Whereas, few other novel miRNAs namely, tpa-miR-N1, tpa-miR-N2, tpa-miR-N6, tpa-miR-N8 and tpa-miR-N9 were found to be very minimal ( 15 copies). The length and Minimum Free Energy (MFE) for these novel pre-miRNAs ranged from 58-82 nts and -20.8 to -44.6 kcal/mol respectively. The (A+U) % of the novel pre-miRNAs was in the range of 37.84% to 60.87% (Table 4).

Sequence and phylogenetic analysis
Sequence and phylogenetic analyses revealed that some of the known miRNAs were expressed in a wide range of insect species and are highly conserved (Table 3, Fig 4A1-4D1 and Fig 4A3-4D3). Mature miRNAs are highly conserved among various species within the Kingdom and are considered to be the evolutionarily conserved regulators of the gene expression [48]. The phylogenetic trees for miR-1000, miR-1175, miR-281 and miR-279 revealed that T. palmi   Discovery of miRNAs from Melon Thrips miRNAs grouped with the closely related species of insects (Fig 4A2-4D2). However, few miR-NAs (miR-1000, miR-2796, miR-965, miR-998, miR-2779, etc.) are highly specific to few species (Table 3). Fig 4A-4D revealed that T. palmi miRNAs are well conserved, particularly in the seed region compared to the homologous miRNAs from other species.

Target Prediction
Targets were predicted for known and novel miRNAs of T. palmi employing miRanda on a scale of 0-7 to indicate the stringency of miRNA-target pairing with the smaller numbers representing higher stringency. ESTs and transcriptome of F. occidentalis were used as a reference for target searches with a cut-off score 140.
Targets for known miRNAs. All 67 known miRNAs were searched for targets against ESTs and transcriptome sequences of F. occidentalis. Out of the known 67 miRNAs, 20 and 40 known miRNAs were found to have targets in ESTs and transcriptome respectively (S2 and S3 Tables). The enrichment analysis (Blast-2-GO) was performed employing gene ontology (GO) terms for genes targeted by miRNAs (Fig 5A and 5B). For those targets in the ESTs, three motifs were over-represented in GO-BP (biological process) category like 'metabolic process' , 'transport' and 'translation' . The GO-MF (molecular function) category was over-represented by the motif 'activity' and 'binding' (Fig 5A). On the other hand, GO-terms enrichment analysis of miRNA targets in the transcriptome yielded motifs for 'transport' , 'metabolic process' and 'oxidation-reduction process' in GO-BP category; while, GO-MF category was over-represented with motifs for 'nucleic acid binding' , 'zinc-ion binding' and ' ATP binding' (Fig 5B). Complete details of the Blast-2-GO analysis have been provided in S4 and S5 Tables.
Targets for novel miRNA. Ten novel miRNAs were searched for their targets in the F. occidentalis transcriptome. A total of 33 miRNA-target pairs were obtained (S6 Table) and   Table). Complete details of the miRNA targets and Blast-2-GO analysis have been provided in S7 Table. The synteny analysis of the T. palmi miRNAs and their targets were performed by employing circos 44 . In brief, the Blast analysis was performed using T. palmi miRNA sequences (known and novel) against F. occidentalis scaffolds (Largest 300). The positions of miRNAs were identified and their targets represented in the Circos plot (Fig 6).

Illumina deep sequencing approach for identification of microRNAs
With the advent of the next generation sequencing technologies, miRNAs have been discovered at an accelerated pace. Presently miRNAs are known from more than 25 insect species, which includes 12 Drosophila species [49]. Among them, the most recent ones are from Plutella xylostella [47], Spodoptera frugiperda [50], etc. Several miRNAs have been reported from various orders of insects such as Diptera, Hymenoptera, Coleoptera, Orthoptera, Lepidoptera, Hemiptera and Homoptera [27], and for the first time, we report the small RNAs from a thysanopteran insect, T. palmi. In this regard, small RNA library was prepared from the pooled samples of different developmental stages of T. palmi and then Illumina (sequencing-by-synthesis) sequencing technology was used to identify miRNAs from the library. The Illumina sequencing approach is one of the high throughput technologies by which miRNAs of any organisms can be identified [51][52][53][54][55][56][57][58][59][60]24]. Size distributions of the high-quality reads varied from 18-28 nts in the library. The peak was at 25 nt which was at par with the previous studies [47,60,61]. According to Bartel (2004), the average miRNA length is 22 nt in animals and study conducted showed that the average length is 23 nt. The unique read distributes of 26-28 nts with a relative lower abundance are common for many small RNA libraries [62][63][64] indicating the presence of piRNAs. Piwi RNAs (piRNAs) are the class of small RNAs mediating chromatin modifications [65] which are derived mainly from retro-transposons and other repetitive elements with high sequence diversity [62,65,66]. Thus, the results indicated that T. palmi genome encodes not only miRNAs but also other statistically significant difference at level p < 0.05 and p < 0.001 respectively for these miRNAs in the larva and adult T. palmi. The error bars indicate standard deviation for three biological replications. (C) Small RNA Northern blot validation. Both conserved (tpa-miR-750, tpa-miR-92b, tpa-miR-2796) and novel (tpa-miR-N3, tpa-miR-N4 and tpa-miR-N10) miRNAs were validated by small RNA northern analysis employing small RNA isolated from pooled T. palmi. doi:10.1371/journal.pone.0163635.g007 Table 6. MicroRNA specific primers employed in the RT-qPCR. Discovery of miRNAs from Melon Thrips small RNAs such as piRNAs (S1 Table) in a lower proportion that might be involved in the trans-generational epigenetic inheritance [67].

Homology-based predictions of miRNAs
The identification of small RNAs (especially miRNAs) based on genomic information has been reported previously in several insects [60,62]. In this regard, we report the identification and characterization of miRNAs from T. palmi based on Illumina small RNA sequencing. We employed F. occidentalis genome sequences as a reference for T. palmi since the complete genome for T. palmi is still not available. However, a large proportion (93.12%) of the T. palmi sRNA sequences could be mapped on to F. occidentalis genome. This higher percentage of mapping was possible because T. palmi & F. occidentalis (reference genome) belong to the same family, Thripidae. Mapping onto a whole genome sequence also helped in elucidating the sample proportion of small ncRNAs such as tRNA, rRNA, snoRNA or snRNA [68]. All of these sequences were annotated by aligning the reads with Rfam database, which indicated the efficiency of deep sequencing in identifying small RNAs. Our results indicated that there is a rich small RNA world evident in T. palmi.
Our study revealed 67 conserved and ten novel miRNAs from T. palmi for the first time. The (A+U) content of the pre-miRNAs should be in the range of 30-70% [69], as those with higher (A+U) content bind more strongly to proteins [20,70]. In this regard, the (A+U) content of the novel pre-miRNAs ranged from 37.84% (tpa-miR-N9) to 60.87% (tpa-miR-N8). In our dataset, the most abundant miRNAs were tpa-miR-281 (from the conserved miRNAs) and tpa-miR-N3 (from the novel miRNAs) with a total of 9560 and 4007 number of reads respectively.
miRNAs are known to be conserved among different species within a kingdom and are evolutionarily conserved regulators of gene expression [20,71]. Our homology and phylogenetic analysis revealed that insect miRNAs are known to be well-conserved, despite considerable diversity in the genome (Fig 4A-4D). In most of the cases, detection of miRNA Ã s is difficult with the available methods as these molecules are liable to degrade soon after being exported to the cytosol [27]. However, in our study during the process of identifying conserved miRNAs, two miRNAs for instance, miR-6489 Ã and miR-6493 Ã that matched to the same precursor sequences with their mismatched complementary mature miRNAs were also detected. The weblogo sequences analysis revealed the likely presence of miR-281 Ã which was not identified in the raw reads of NGS. The presence of miR-281 Ã was identified by BLASTN option in miR-Base and further confirmed by stem-loop RT PCR (Table 5, Fig 4B3). The absence of miR-281 Ã could be due to the faster degradation as compared to miR-281.

Possible roles of T. palmi miRNAs
Although thousands of small RNAs have been discovered in the recent past, [39,40,47,50,72] the primary challenge is to fully identify the spatiotemporally expressed microRNAs and to determine their individual functions. The majority of the microRNAs have been identified through either computational prediction or cloning and sequencing [73]. In this study, we employed Illumina next generation sequencing approach to identify miRNAs from T. palmi. Currently, there are several mature and precursor microRNAs deposited in the miRBase [39]. In this connection, we identified a total of 77 miRNAs from T. palmi using high throughput sequencing. The current study is the first report of miRNA profiling from a Thrips species employing deep sequencing approach. This approach is far superior to the other approaches of miRNA identification, as it can discover novel microRNAs [74].
The analysis of the expression value (read numbers) revealed that the highest expression was for miR-281 (9560 reads). Recent studies have proved that microRNA-281 regulate the expression of ecdysone receptor (EcR) isoform B, in the silkworm, Bombyx mori [75]. Thus, miR-281 may be involved in development and metamorphosis of T. palmi by regulating the genes involved in the ecdysone cascade. The second highest expression was for miR-750 with an expression value of 7869. RNAi studies proved that the putative JH receptor ultraspiracle (USP) [76,77] is a likely target of miR-750. Thus, it indicates that miR-750 may be involved in hormone signaling, immunity and stress response by regulating the vitellogenin (Vg) gene in T. palmi [78].
Another interesting microRNA obtained in the current study was miR-92 with an expression value of 1687. Previous studies had shown that miR-92 regulates Mef2, the key transcription factor for muscle development and differentiation in Drosophila [79]. An insect-specific microRNA, miR-2796 was identified in T. palmi with an expression value of 479. miR-2796 was found to be the most abundant microRNA in honey bee brain [80]. Additionally, miR-2796 bound to the coding sequence (CDs) of PLC-epsilon gene in Apis and Tribolium, affecting the mRNA stability by splicing rather than the normal canonical translational repression [81]. However, interestingly both miR-2796 and PLC-epsilon gene were missing from the genus Drosophila, even though it was found in other dipterans [80].
Our analyses revealed miR-993 were identified only in invertebrates (Table 3). miR-993 belongs to the miR-100/10 family, and both miR-993 and miR-10 are derived from the ancient miR-100 through duplication and arm-switching [60,82]. miR100/10 family members could regulate the expression of relevant Hox-genes, thus may play a crucial role in insect development [83]. Rest of the insect-specific miRNAs identified in T. palmi may play an important role in insect-specific features, such as metamorphosis, parthenogenesis and biogenesis of pheromones [84]. Whereas, the other invertebrate and vertebrate-specific miRNAs (Table 3) identified from T. palmi required special attention, as their nonexistence in other species of insects could be due to the absence of genomic information for most of those insects [60]. These specific microRNAs may be involved in some special biological processes that distinguish thysanopteran insects from others.
Target prediction is crucial to understand the biological role of a particular miRNA. Unlike their plant counterparts, the imperfect complementarity of animal miRNAs to their target mRNA sequences makes it more difficult to judge the accuracy of prediction. MicroRNAs can bring about mRNA cleavage or translational repression of target mRNAs by binding to 3' UTRs, 5' UTRs and even to coding regions [87]. However, animal miRNAs primarily target the 3' UTRs; and therefore, we limited our target search to (i) expressed sequence tags (ESTs) and (ii) transcriptomic sequences of F. occidentalis. The predicted targets were annotated against GO database and the targeted genes included transcription factors, signal transduction, hormone pathways, molting and even metabolism. Therefore, all these conserved and novel miRNAs identified from T. palmi could play a vital role in diverse biological processes, thus undoubtedly participating in the regulation of thrips growth and development.
In summary, results from this study add to our growing pool of miRNA database and is the first report on such analysis in a thysanopteran insect, T. palmi. Deep sequencing of small RNAs has facilitated the identification of miRNAs from T. palmi. Sixty-seven conserved and ten novel miRNAs that were identified with high confidence and sufficient evidence are the contributions from our study. Most of the T. palmi miRNAs were homologous to insects as compared to the vertebrates. Sequence and phylogenetic analyses revealed that most of the T. palmi miRNAs are highly conserved in various species, making miRNAs, a hallmark of evolutionarily conserved regulators of gene expression. To harmonize the data, and to provide more useful biological insights, we also carried out in silico analysis for identifying potential targets for these miRNAs. Unlike the plant counterparts, the imperfect complementarity of metazoan miRNAs to the target has been found to be sufficient to promote the RNA silencing, as in the case of Drosophila and Bactrocera [88]. Our results indicated that the list of putative mRNA targets was very extensive (S2, S3 and S6 Tables), even with stringent parameters applied to miRanda. Our results suggested that most of the putative target genes for T. palmi miRNAs were associated with several KEGG pathways like metabolic process, transport, translation, signal pathways and oxidative phosphorylation. However, further wet lab experiments are still required for the validation of these targets in understanding the biology of this insect. Expression levels of few miRNAs were also validated by both qRT-PCR and Northern analysis.
Several miRNAs were identified and characterized from animals and plants and among them, very few were further explored for various applications by disrupting specific pathways targeted by these miRNAs. This can be achieved by employing the artificial microRNAs (amiRs) [89]. Recent studies successfully demonstrated the use of amiRNAs for targeting the reporter and the endogenous genes in animals and plants. Identification and expression of a few essential insect-specific gene(s) in plants, can target and degrade an invading insect's genes, consequently confer insect resistance [90]. Results from our study indicated few miR-NAs have been predicted to be involved in the adult development process, which can be further utilized in gene functional studies through RNAi-based approach or in developing miRNA mimics both for feeding and in planta expression [29,88,91] as novel pest management strategies based on gene silencing and insect transgenesis.
Supporting Information S1