Genome‐wide identification, expression profiling, and target gene analysis of microRNAs in the Onion thrips, Thrips tabaci Lindeman (Thysanoptera: Thripidae), vectors of tospoviruses (Bunyaviridae)

Abstract Thrips tabaci Lindeman is an important polyphagous insect pest species estimated to cause losses of more than U.S. $1 billion worldwide annually. Chemical insecticides are of limited use in the management of T. tabaci due to the thigmokinetic behavior and development of resistance to insecticides. There is an urgent need to find alternative management strategies. Small noncoding RNAs (sncRNAs) especially microRNAs (miRNAs) hold great promise as key regulators of gene expression in a wide range of organisms. MiRNAs are a group of endogenously originated sncRNA known to regulate gene expression in animals, plants, and protozoans. In this study, we explored these RNAs in T. tabaci using deep sequencing to provide a basis for future studies of their biological and physiological roles in governing gene expression. Apart from snoRNAs and piRNAs, our study identified nine novel and 130 known miRNAs from T. tabaci. Functional classification of the targets for these miRNAs predicted that majority are involved in regulating transcription, translation, signal transduction and genetic information processing. The higher expression of few miRNAs (such as tta‐miR‐281, tta‐miR‐184, tta‐miR‐3533, tta‐miR‐N1, tta‐miR‐N7, and tta‐miR‐N9) in T. tabaci pupal and adult stages reflected their possible role in larval and adult development, metamorphosis, parthenogenesis, and reproduction. This is the first exploration of the miRNAome in T. tabaci, which not only provides insights into their possible role in insect metamorphosis, growth, and development but also offer an important resource for future pest management strategies.

Thrips tabaci is also a vector of two viral pathogens, Iris yellow spot virus (IYSV) (Srinivasan et al., 2012) and Tomato spotted wilt virus (TSWV) (Pittman, 1927) causing significant disease around the world (German, Ullman, & Moyer, 1992). Thrips tabaci is estimated to cause more than U.S. $1 billion in crop losses annually worldwide. To date, chemical insecticides have been widely used for the management of T. tabaci, but due to its thigmokinetic behavior and frequent development of insecticide resistance, they have had little use. Therefore, the design of novel insecticides, resistance breeding strategies, an in-depth understanding of genes and gene regulation is necessary for targeting important developmental factors/processes for effective management of this insect. MiRNA analysis is an effective tool to understand gene regulation and expression in both insect and host plant.
So far, the miRNAome for insects is far behind nematodes, plants, and mammals (Kakumani et al., 2015). MiRNAs are reported from about 25 species of insects belonging to various orders (Stark et al., 2007;Wu et al., 2013). No information is available on T. tabaci miRNA content and function. Our study reports the detailed profile of miRNAs from T. tabaci. Further analysis identified putative target genes for these miRNAs, which will shed more light on the identification of highly specific miRNAs for thysanopteran pest management in the near future.

| Insect culture and RNA isolation
Thrips tabaci cultures were maintained on Phaseolus vulgaris in controlled laboratory conditions at 25°C (DeGraaf & Wood, 2009) with an 8 hr:16 hr light:dark cycle. Total RNA was isolated from whole-body homogenates of sample mix, containing a total of 50 mg of different life stages viz. eggs, larvae, pupae, and adults of T. tabaci using TRIzol reagent (Invitrogen, Carlsbad, CA, USA).

| Sample preparation and Illumina sequencing
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.
Homology analysis was performed with conserved miRNAs of T. tabaci with the miRNAs of other organisms from the miRBase database (Release 21.0; Griffiths-Jones, Saini, van Dongen, & Enright, 2008). BLASTn embedded in the miRBase database was used to compare the T. tabaci miRNAs with other species, with an E-value of .01 to find out more miRNA homologs. The naming of the miRNAs in this study has been performed according to Griffiths-Jones, Grocock, van Dongen, Bateman, & Enright, 2006. As these miRNAs were predicted from T. tabaci, the prefix for all miRNAs was fixed as "tta." The rest of the naming convention criteria were in accordance with miRBase (Griffiths-Jones et al., 2006).

| Target prediction
Targets for identified miRNAs were predicted employing the miRanda program (Enright et al., 2004), against the expressed sequence tags (ESTs) and transcriptome (NCBI Accession: PRJNA203209) database of Frankliniella occidentalis. An alignment score (Smith & Waterman, 1981)  T A B L E 3 Expression value of known miRNAs in Thrips tabaci. The first column represents miRNA family; the second column represents the number of reads annotated on the particular miRNA family; the third column represents length of the mature miRNA sequences; the fourth column represents the name of the miRNAs in T. tabaci; the fifth column represents the miRNA sequence; the sixth column represents homologous species of organism where it has the highest similarity  were assigned (using Blast-2-GO) based on the annotation. The circos plot was generated using Circos (Krzywinski et al., 2009) to visualize the interaction between miRNAs and their targets.

| Validation of Thrips tabaci miRNAs using Stemloop RT-PCR
We were able to validate six conserved and four novel microRNAs employing Stem-loop RT-PCR primers designed based on previous reports (Chen et al., 2005).  (Bustin et al., 2009). qRT-PCR assays were performed in triplicates for three independent biological replicates, and the relative gene expression data were analyzed using 2 −ΔΔC T method (Livak & Schmittgen, 2001). U6 snRNAs was used as an internal control gene for normalization. The values of these three independent experiments were statistically analyzed using one-way ANOVA to calculate the statistical significance.

| Illumina sequencing of Thrips tabaci small RNAs
The small RNA library prepared for deep sequencing resulted in a total of 13,192,454 raw reads ( both known and novel miRNAs. Size distributions of the trimmed highquality reads were varied from 18 to 26 nts with a peak at the 23 nts ( Figure 2). A small portion of our library consisted of read length of around 26-28 nts, which could be putative piwi-interacting RNAs (pi-RNAs) from T. tabaci as the homology search against the piRNABank database revealed that some of these were similar to previously reported piRNAs (Table 2).

| Identification of known miRNAs from Thrips tabaci
Our analyses on the trimmed high-quality reads resulted in a total of 130 conserved miRNAs representing 55 different miRNA families (Table 3).

| Identification of novel miRNAs from Thrips tabaci
Miranalyzer pipeline identified a total of nine novel miRNAs from T. tabaci for the first time (Table 5), with their predicted precursor secondary structures (Figure 3). The complete details of the mature miRNAs and their corresponding pre-miRNAs have been given in  Figure 3).

| The presence of miRNA star strands
It is very difficult to identify the star strand (miRNA*) sequences from the library, as it will be degraded soon after being exported to the cytosol. However, our results revealed that ten also indicated the presence of miRNA* sequences in four of our novel miRNAs such as tta-miR-N6, tta-miR-N7, tta-miR-N8, and tta-miR-N9, although the abundance was low ( Table 5). The complete characteristic features of these miRNA* sequences and their corresponding pre-miRNA*s have been given in Tables 5 and 6.

| Identification of plant miRNA family in Thrips tabaci sRNA library
Interestingly, this study has identified mir-9774 (Expression value 6), a plant microRNA family in our T. tabaci sRNA library (Table 3).

| Phylogenetic analysis of Thrips tabaci miRNAs
Phylogenetic analyses revealed that most of the known miRNAs are highly conserved (

| Identification of targets for Thrips tabaci miRNAs
Targets were predicted for known and novel miRNAs of T. tabaci employing miRanda with 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 from Thrips tabaci
One hundred and thirty known miRNAs were searched for targets against ESTs and transcriptome sequences of F. occidentalis. A total of 218 and 1,025 targets were obtained from ESTs and transcriptome, respectively (Tables S1 and S2). The Blast-2-GO enrichment analysis was performed employing gene ontology (GO) terms for genes targeted by these miRNAs (Figure 5a,b). For those targets in the ESTs, three motifs were over-represented in GO-BP (biological process) category viz.
"metabolic process," "transport," and "catabolic process." The GO-MF (molecular function) category was over-represented by the motif "oxidoreductase activity" and "catalytic activity" (Figure 5a). On the other hand, GO terms enrichment analysis of miRNA targets in the transcriptome yielded motifs for "transport," "signal transduction," and "metabolic process" in GO-BP category; while, GO-MF category was over-represented with motifs for "ATP binding," "transferase activity," and "binding" (Figure 5b). Complete details of the Blast-2-GO analysis were provided in Tables S3 and S4.

| Targets for novel miRNAs from Thrips tabaci
Novel miRNAs were searched for their targets in the F. occidentalis transcriptome. A total of 65 miRNA-target pairs were obtained (Table   S5), and further Blast-2-GO analysis indicated the over-representation of "Transport" and "ATP binding" as GO-BP and GO-MF category, respectively ( Figure 6 and Table S6).

| Synteny analysis using Circos
The synteny analysis of the T. tabaci miRNAs and their targets were performed by employing circos (Krzywinski et al., 2009). In brief, the

| Validation of Thrips tabaci microRNAs
This study revealed 130 known and nine novel miRNAs from T. tabaci. However, further validation of these miRNAs was performed by (1) stem-loop endpoint reverse transcriptase PCR (RT-PCR) and (2) real-time quantitative reverse transcriptase PCR (RT-qPCR). Using stem-loop endpoint RT-PCR, we have validated six conserved viz.
F I G U R E 4 (a-d): 1. Homology in the seed region of the Thrips tabaci miRNAs (a-d are for mir-8, mir-14, mir-276, and mir-281, respectively) with respect to its counterpart from other insect species. The first three letters of each miRNAs indicating the name of the species (e.g.,: dya-Drosophila yakuba).

| DISCUSSION
The onion thrips, Thrips tabaci, is an important pest species and a tospovirus vector causing significant negative impacts on yield and quality of various economically important crops (German et al., 1992).
Although microRNAs are key gene regulators and are involved in many biological processes, including growth and development, no previous study has been conducted on the identification and validation of miRNAs in T. tabaci. MicroRNAs are known from more than 25 insect species, (Stark et al., 2007). Several miRNAs have been reported from various orders of insects such as Diptera, Hymenoptera, Coleoptera, Orthoptera, Lepidoptera, Hemiptera, Homoptera (Wu et al., 2013), and Thysanoptera (Rebijith, Asokan, Hande, & Krishna Kumar, 2016). This study reports the complete miRNA profile from onion thrips, Thrips tabaci. A small RNA library was prepared from the pooled samples of different developmental stages of T. tabaci and the high-throughput Illumina deep-sequencing technology (Avesson et al., 2012;Burnside et al., 2008;Ge et al., 2013;Koh et al., 2010) was used to identify miRNAs from the prepared library.
We used the F. occidentalis genome sequence as a reference for T. tabaci, as the complete genome T. tabaci is still not available in the database. The higher percentage of mapping (91%) was possible only because both these insects belong to the same family, Thripidae.
Employing this approach, our study revealed 130 conserved and nine novel miRNAs from T. tabaci. The size distributions of the highquality reads were varied from 18 to 28 nts in our library and the peak was at the 25 nt, which was on par with previous studies Liang, Feng, Zhou, & Gao, 2013;Sattar et al., 2012). Our study indicated the unique read distributes of 26-28 nts with a relative lower abundance, which is common in many small RNA libraries (Chang et al., 2016;Jagadeeswaran et al., 2010;Surridge et al., 2011;Zhang et al., 2013), indicating the presence of piRNAs. Piwi RNAs (piRNAs) are the class of small RNAs mediating chromatin modifications (Ross, Weiner, & Lin, 2014) which are derived mainly from retrotransposons and other repetitive elements with high sequence diversity (Ross et al., 2014;Siomi, Sato, Pezic, & Aravin, 2011;Zhang et al., 2013). Thus, our results indicated that T. tabaci genome not only harbors miRNAs but also other small RNAs such as piRNAs that might be involved in the transgenerational epigenetic inheritance (Weick & Miska, 2014).
MiRNAs are evolutionarily conserved regulators of gene expression (Rebijith et al., 2014;Zhang et al., 2009), and few can even act as markers in defining the evolutionary relationship in a wide range of insect species (Kakumani et al., 2015). Our homology and phylogeny analysis revealed that insect miRNAs are well-conserved, despite considerable diversity in the genome (Figure 4a-d). MiRNA*s are not easily detectable as it degrades soon after being exported to the cytosol (Wu et al., 2013). However, our results indicated the presence of several miRNA*s (Tables 5 and 6) that matched to the same precursor sequences with their mismatched complementary mature miRNAs.
We identified the presence of a plant-specific miRNA family, mir-9774 in the T. tabaci sRNA library, and the same has been recently reported from Triticum aestivum L. and Brachypodium distachyon (L.) Beauv (Wei et al., 2009). Previous miRNA studies on cotton/melon aphid, A. gossypii also reported six plant miRNA family (Sattar et al., 2012). They also showed that such microRNAs were transformed into the aphid tissues (especially in gut contents) during the phloem sap ingestion. However, none of those six have been identified in our sRNA library.
Our results showed that the highest expression is for tta- is quite possible that miR-8 may play a key role in the reproductive processes of T. tabaci. An insect-specific miR-14 was identified in T. tabaci with an expression value of 12,453 and studies on lepidopteran insects showed the antiapoptotic role of miR-14 (Kumarswamy & Chandna, 2010). The rest of the species-specific miRNAs identified in T. tabaci might play important role in insectspecific features, such as metamorphosis, parthenogenesis, and biogenesis of pheromones (Zhang et al., 2007). Whereas, the other invertebrate-and vertebrate-specific miRNAs (Table 3) identified from T. tabaci required special attention, as their nonexistence in other species of insects could be due to the absence of complete genomic information for most of those insects .
The expression profile of miRNA varies spatiotemporally among different developmental stages (Li, Cassidy, Reinke, Fischboeck, & Carthew, 2009;Xu, Zhou, Wang, Auersperg, & Peng, 2006), and the developmental expression profiles (larval, pupal and adult stage) of scripts that can bring about mRNA cleavage, mRNA decay or translational repression of target mRNAs by binding to 3′ UTRs, 5′ UTRs, and even to coding regions (Lytle, Yario, & Steitz, 2007). Thus, it is important to identify the gene targets and thereby we can understand the biological role of a particular miRNA. As miRNA targets have been identified using the (1) expressed sequence tags (ESTs) and (2) transcriptomic sequences of F. occidentalis. The GO annotations for the predicted targets were classified as potential biological process, cellular component, and molecular function. The putative targeted genes included signal transduction pathways, transcription factors, reproduction, embryo development, insect molting, immune response, and even metabolism. Overall, the results from our study indicated that these conserved and novel miRNAs identified from T. tabaci might play crucial regulatory role in the regulation of thrips growth and development.
F I G U R E 7 Map of the Western Flower Thrips, F. occidentalis scaffolds linking T. tabaci miRNAs and their putative targets prepared using Circos (Krzywinski et al., 2009). The outer circle represents the highlights of nine novel miRNA represented in blue and 34 known miRNA represented in red color. The inner circle marks each scaffold in a different color. The blue lines in the center of the figure connect a known miRNAs with its target that are represented across 173 scaffolds of F. occidentalis genome. Whereas, the orange lines in the center represent the interaction of novel miRNA with its target positions

| CONCLUSIONS
In summary, the result from our study add to the pool of miRNA databases and is the first report of small RNAs from T. tabaci, a nonmodel insect lacking genome information. One hundred and thirty conserved and nine novel miRNAs were identified with high confidence and sufficient evidence is the major contribution of our study. Sequence analyses revealed that most of the T. tabaci 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 have also carried out in silico analysis of identifying potential targets for these miRNAs. Our results indicated that the list of putative mRNA targets was very extensive and most of the putative target genes for T. tabaci miRNAs were associated with several KEGG pathways such as metabolic process, transport, translation, signal pathways, and oxidative phosphorylation. However, further experiments are required for the validation of these targets. Expression levels of T. tabaci miRNAs were validated by RT-qPCR, and the results indicated few of these miRNAs have been predicted 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 (Agrawal, Sachdev, Rodrigues, Sowjanya Sree, & Bhatnagar, 2013;Jayachandran, Hussain, & Asgari, 2013;Nandety et al., 2015) as novel pest management strategies based on gene silencing and insect transgenesis.

ACKNOWLEDGMENTS
We thank Sonal Dsouza, Communication and Programme Assistant; Bioversity International for language editing that greatly improved the manuscript. Our sincere thanks are due to The Director, IIHR, Bangalore for providing necessary facilities.

DATA AVAILABILITY
All relevant data are within the paper and its Supporting Information files. The small RNA Sequence data has been submitted to NCBI under the BioSample project 'PRJNA350618'; BioSample Accession: 'SAMN05943039'.

CONFLICT OF INTEREST
The authors have declared that no competing interests exist. F I G U R E 8 (a) Stem-loop RT-PCR analyses of six conserved and four novel miRNAs from Thrips tabaci. The products were resolved on 3% agarose gel in 1X TBE stained with ethidium bromide and HyperLadder ™ 25 bp (Bioline, USA) used as marker. (b) Stem-loop RT-qPCR analysis of spatiotemporally expressed T. tabaci miRNAs in larva, pupa and adults. "*" and "**" means a statistically significant difference at level p < .05 and p < .001, respectively, for these miRNAs in the larva, pupae, and adult T. tabaci. The error bars indicate standard deviation for three biological replications