Characterization and Comparative Profiling of miRNAs in Invasive Bemisia tabaci (Gennadius) B and Q

Background MicroRNAs (miRNAs) are small, conserved, non-coding RNAs that post-transcriptionally regulate gene expression. Bemisia tabaci (Gennadius) B and Q are two invasive and dominant whiteflies, and B. tabaci Q has been displacing B in China. Differences in biological traits (fecundity, host range, resistance to insecticides, etc.) as affected by miRNAs might be involved in the displacement. In this study, we performed high-throughput sequencing to identify miRNAs in B. tabaci B and Q. Results We identified 170 conserved miRNAs and 15 novel candidates, and found significant differences in the expression of miRNAs between B. tabaci B and Q. Conclusion Expression levels of miRNAs differ in B. tabaci B vs. Q. Additional research is needed to determine whether these differences are related to differences in the biology of B. tabaci B and Q, and whether these differences help explain why B. tabaci Q is displacing B in China.

Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae) is an important agricultural pest worldwide that attacks more than 600 plant species including food, fiber, and ornamental plants under field and greenhouse conditions [32,33]. By phloem feeding, contaminating leaves and fruits with honeydew, and transmitting more than 110 kinds of plant viruses, adult and immature instars of B. tabaci cause billions of dollars of annual loss worldwide [34,35]. In particular, B. tabaci often causes outbreaks of plantpathogenic viruses [36,37].
Bemisia tabaci is currently regarded as a cryptic species complex that contains at least 24 morphologically indistinguishable species [38]. Two members of this complex, Middle East-Asia Minor 1 (commonly known as B. tabaci biotype B, herein referred to as B. tabaci B) and Mediterranean (commonly known as B. tabaci biotype Q, herein referred to as B. tabaci Q), have now spread well beyond their home ranges as a consequence of trade in ornamental plants [38]. In China, B. tabaci B and Q are the main whiteflies in agricultural areas and have caused severe damage to many crops [39][40][41]. We previously reported that the ratio of abundances of B. tabaci Q to B. tabaci B has been increasing and that Q is displacing B on cotton, eggplant, and other plants in Shandong Province of China [42,43]. Teng et al. [40] also found that B. tabaci Q has become dominant across China and suggested that B. tabaci Q has been displacing non-Q whiteflies in many regions such as Shanxi, Henan, Hubei, Jiangsu, Zhejiang, Hunan, and Hainan provinces.
As noted earlier, many biological differences between B. tabaci B and Q have been documented. B. tabaci Q had significantly greater reproductive parameters than B in winter weeds [44], had shorter developmental times than B on sweet pepper at constant temperatures [45], had better survival than B under low as well as high temperature conditions [46], and had greater resistance to neonicotinoides and pyriproxyfen insecticides than B [47][48][49]. All of these would seem to at least partially explain the displacement of B by Q [50]. The genetic differences between B and Q have also been studied recently, and several genes involved in metabolism and insecticide resistance were considered as possibly contributing to the divergence of the two whitefly species [51]. Because miRNAs are now recognized as critical regulators of gene expression and animal development, the identification and comparison of miRNAs in B. tabaci B and Q could provide new information about the biological differences in these biotypes. The differences between miRNAs of B and Q might also explain the biological differences because miRNAs play an important role in a wide range of cellular and developmental process including cell proliferation [52], cell differentiation [53,54], the cell cycle [55,56], metabolism [57], developmental timing [58], reproduction [8], apoptosis [59], and others [11].
In recent years, high-throughput sequencing technology has been widely used to analyze the characteristics of miRNA within the organisms [60][61][62]. For example, by high-throughput sequencing of miRNA, Marco et al. [60] characterized 203 miRNAs from the red flour beetle Tribolium castaneum; Liu et al. [61] analyzed areas of skin where the cashmere grows in anagen and found that the miRNAs that were coexpressed in goat and sheep were located in the same region of the respective chromosomes and may play an essential role in skin and follicle development; Shao et al. [62] analyzed Arabidopsis (Arabidopsis thaliana) and rice (Oryza sativa) and found that the accumulation levels of several miRNA*s could be much higher than those of their miRNA partners in certain organs, mutants and/or AGOassociated silencing complexes of both Arabidopsis and rice.
In this study, we performed high-throughput sequencing to identify miRNAs in B. tabaci B and Q, and identified 170 conserved miRNAs and 15 novel candidates. We compared the expression of miRNAs in B and Q to identify differentially expressed miRNAs. The results indicate significant differences in the expression of miRNAs between B. tabaci B and Q.

Whitefly Colony and RNA Extraction
Bemisia tabaci B and Q colonies were reared on cotton leaves in growth chambers at 2661uC and with a 16/8 h light/dark photoperiod. Adult whiteflies (B. tabaci B and Q) were collected and homogenized in Trizol agent RNAiso Plus (TaKaRa, Dalian, China). Total RNA was extracted from B and Q according to the manufacturer's instructions and was quantified with an Agilent 2100 Bioanalyzer.

Confirmation of B. tabaci B and Q
The identities of B. tabaci B and Q were confirmed based on the cleaved amplified polymorphic sequences (CAPS) of mtCOI (mitochondrial cytochrome oxidase I) with the restriction endonucleases VspI and StuI [43,63]. Genomic DNA was extracted from individual adult whiteflies according to the lysis method of Frohlich et al. [64]. The mtCOI fragments were amplified using primers C1-J-2195 (59-TTGATTTTTTGGTCATCCAGAAGT-39) and R-BQ-2819 (59-CTGAATATCGRCGAGGCATTCC -39) [63]. The 20 mL PCR reaction mixture contained 2 mL of 106reaction buffer supplemented with 1.5 mM MgCl 2 , 0.2 mM of each primer, 0.2 mM of each dNTP, 1 unit of Taq DNA polymerase, and 2 mL of each template cDNA. Cycling conditions were as follows: 5 min at 94uC; 35 cycles of 1 min at 94uC, 1 min at 52uC, and 1 min at 72uC; and finally 10 min at 72uC. PCR products were electrophoresed and visualized by ethidium bromide staining. The mtCOI fragment (approximately 620 bp) was first cleaved by VspI [65], and then the uncut fragment was cleaved by StuI [66]. Specimens whose mtCOI fragments were cut by VspI were identified as B. tabaci Q, whereas specimens whose mtCOI fragments were cut by StuI were identified as B. tabaci B [43,63].

Small RNA Library Preparation and High-throughput Sequencing
For HiSeq deep sequencing, the small RNA samples were prepared as described previously [29]. In brief, RNA fragments with fewer than 40 nt were isolated from total RNA on a 15% Novex TBE-urea PAGE gel. Then, a 59 adaptor (Illumina, San Diego, CA, USA) was ligated to the purified small RNAs, and the ligation products were purified on a 15% Novex TBE-urea PAGE gel. The 59 ligation products were then ligated to a 39 adaptor (Illumina), and products with 59 and 39 adaptors were sizefractionated on a 10% Novex TBE-urea PAGE gel. Subsequently, small RNAs ligated with adaptors were reverse transcribed and then subjected to PCR amplification. The amplification products were purified on a 6% Novex TBE PAGE gel. The purified DNA fragments were used for clustering and sequencing by HiSeq highthroughput sequencing technology at the Beijing Genomics Institute, Shenzhen.

Discovery of Conserved miRNAs
The tags under 40 nt sequence from HiSeq sequencing were first subjected to data cleaning, which included removal of the low quality tags and several kinds of contaminants. The distribution of the lengths of the clean tags was then summarized, and the clean tags were assigned to two groups including the summary of unique tags and total tags. The clean tags were annotated into different categories to discard rRNAs, tRNAs, snRNAs, and snoRNAs using Rfam database (version 10.1). Because there was no information concerning miRNAs of B. tabaci in the miRBase v17.0, the remaining small RNA tags were aligned to the miRNA precursors/mature miRNAs of all animals in the miRBase v17.0 [67,68]. Sequences in our libraries that were identical to or related to (having four or fewer nucleotide substitutions) sequences from Drosophila melanogaster or other insects (A. aegypti, A. mellifera, B. mori, and T. castaneum) were identified as conserved miRNAs.

Prediction of Novel miRNA Candidates
The characteristic hairpin structure of miRNA precursors can be used to predict novel miRNA candidates. Because there were no completed genome sequences, 27,288 nucleotide sequences of Bemisia tabaci obtained from NCBI (published by Zhejiang University) were used as a reference for novel miRNA prediction. The prediction software Mireap was used to predict novel miRNA candidates by exploring the secondary structure, the Dicer cleavage site, and the minimum free energy of the unannotated small RNA tags that could be mapped to the genome. The rules used to identify novel miRNA candidates were based on those suggested by Allen et al. [69] and Schwab et al.

Comparing the Expression of miRNAs between B. tabaci B and Q
We compared the expression of miRNAs between B. tabaci B and Q to identify differentially expressed miRNAs. The expression of miRNAs in the two libraries was visualized on a scatter plot in which expression of B miRNAs was plotted against expression of Q miRNAs after expression levels were normalized and then transformed into fold-change values (see below). The threshold of a fold change more than 2 was considered significant difference. The procedure had two parts. First, the expression of miRNA in the two libraries (Q as control and B as treatment) was normalized    to transcripts per million (TPM). If the normalized expression of a miRNA was zero, it was modified to 0.01 to enable calculation. If the normalized expression of a miRNA was less than 1 in both B and Q libraries, it was ignored to compare for low expression. The normalization formula was: Normalized expression = Actual miRNA count/Total count of clean reads*1000000. Second, the normalized data were used to calculate fold-change values and P-values, and a scatter plot of the fold-change values was generated. Fold-change was calculated as: Fold-change = log 2 (B/Q). The P-value was calculated as x represents Q; y represents B; N1 represents the normalized expression of a miRNA in Q library; N2 represents the normalized expression of the same miRNA in B library.

B. tabaci has a Complex Population of Small RNAs
HiSeq high-throughput sequencing technology was used to identify miRNAs in B. tabaci B and Q. Two libraries of small RNAs were constructed, one from B and the other from Q. We obtained 17,953,732 reads from the B library, and 16,448,832 reads from the Q library. Low-quality sequences and those shorter than 18 nt were removed, leaving 17,451,513 reads (2,871,240 unique sequences) in the B library and 15,977,474 reads (3,041,144 unique sequences) in the Q library. The distribution of sequence lengths indicated that both libraries were enriched with small RNAs of 21-23 nt (42.8% and 32.6% of all reads in B and Q libraries, respectively) ( Fig. 1), which is considered the standard size of miRNAs. Another type of RNA sequence found in both libraries was 28-30 nt long, corresponding to pi-RNA-like sequences, and represented 28.4% and 49.4% of the reads in B and Q libraries, respectively. In these libraries, sequence length was limited to a maximum of 30 nt.
Subsequent sequence analysis eliminated reads corresponding to rRNAs, tRNAs, snRNAs, and snoRNAs. Another two large fractions of reads were derived from unannotated genome sites

Identification of Conserved miRNAs
To identify conserved miRNAs in B. tabaci, we searched all small RNA sequences against the currently miRNAs of all animals in miRBase v17.0 using BLASTn [67,68]. In total, 1,504 miRNAs were found in the B library, and 1,182 miRNAs were found in the Q library. Sequences in our libraries that were identical to or related to (having four or fewer nucleotide substitutions) sequences from D. melanogaster or other insects (A. aegypti, A. mellifera, B. mori, and T. castaneum) were identified as conserved miRNAs. After BLASTn searches and further sequence analysis, a total of 170 conserved miRNAs identified from the miRNAs were found in both B and Q libraries (Table 1). In these conserved miRNAs, there were nine miRNA families containing five or more miRNAs ( Table 1). The identified miRNA families are conserved in a variety of animal species. For example, let-7 [71], miR-9 [72,73], miR-10 [74], miR-133 [75,76], and miR-263 [6] have been found in 100, 63, 46, 89, and 74 animal species, respectively.

Identification of Candidate Novel miRNA Candidates
In addition to the identification of conserved miRNAs, we identified 15 potential novel miRNA candidates in both B and Q libraries ( Table 2). The length of the 15 predicted novel miRNA candidates ranged from 21 to 24 nt. The free energy of folding for these hairpin structures ranged from 232.22 kcal/mol to 218.6 kcal/mol. The read number for each novel miRNA was lower than that of the conserved miRNAs, which was consistent with previous studies. To investigate the conservation of these novel miRNA candidates in a wide range of animal species, we used these miRNAs as query sequences to perform BLASTn searches against all nucleotide sequences in miRBase v17.0 databases. No homologs were found in any animal species, suggesting that these newly identified miRNAs are all whiteflyspecific.

Comparing the Expression of miRNAs between B and Q
To explore their difference in miRNAs, we compared the expression of miRNAs between B and Q ( Fig. 3 and Table S1). The expression of 342 miRNAs was higher in the Q library than in the B library, and 198 miRNAs were found only in the Q library. The expression of 474 miRNAs was lower in the Q library than in the B library, and 303 miRNAs were found only in the B library. About 579 miRNAs were found in both B and Q libraries; among these, the expression of 144 miRNAs was higher in the Q library than in the B library, and the expression of 171 miRNAs was lower in the Q library than in the B library.
Discussion miRNAs post-transcriptionally regulate gene expression by targeting the 39 untranslated region of specific messenger RNAs and causing mRNA degradation or translational repression [77,78]. In this study, the expression of the most conserved miRNAs was 1.1 to 2.5 times greater in B. tabaci B than in Q (Table 1). Because miRNAs could recognize the target mRNAs based on sequence complementarity and cause mRNA degradation, the expression of some functional proteins should be lower in B. tabaci B than in Q, which possibly contributing to the biological differences between B and Q.
We found substantial differences between B. tabaci B and Q in the expression of miRNAs. For example, miR-139*, miR-1468, miR-4496, miR-1566, and miR-2993 were found only in B. tabaci  [47][48][49]. Accumulating evidence indicates an important role of miRNAs in drug resistance, and miRNA expression profiling is correlated with the development of drug resistance [79][80][81]. Although many studies have reported the involvement of miRNAs in drug resistance, few of these have concerned insects. We therefore analyzed the miRNAs possibly related to drug resistance on the basis of the studies in humans. In the present study, the normalized expression of miR-638 was 19,471 times higher in B. tabaci B than in Q (Table S1), suggesting a much lower level of drug resistance in B than in Q. Haenisch et al. [82] found that miR-379 increased (maximally 4.1061.33-fold) in HepG2 cells after 48 h of treatment with 5 mM rifampicin. In our study, the normalized expression of miR-379 was 58 times higher in B. tabaci Q than in B (Table S1), again indicating a much greater drug resistance. Additional research is needed to determine whether miRNAs are involved in insecticide resistance.
This miRNA of miR-146b, which is highly homologous to miR-146a, was much more abundant in B. tabaci Q than B, and miR-146c was found only in Q (Table S1). The normalized expression of miR-215 was lower in B. tabaci Q than in B (Table S1), indicating lower apoptosis in Q. Apoptosis is considered a vital component of various processes including normal cell turnover, embryonic development, and chemical-induced cell death [83]. To date, several miRNAs have been identified that regulate the complex networks of apoptotic pathways [84]. Experimental evidence in human has demonstrated that miR-146a modulates activation-induced cell death (AICD), acting as an anti-apoptotic factor, and that was associated death domain (FADD) is a target of miR-146a [85]. Differences in the expression of these miRNAs might also explain difference in survival among B. tabaci biotypes. The tumor suppressor p53 acts as a major defense against cancer and can elicit both apoptotic death and cell cycle arrest [86]. miR-215 was identified as a p53-regulated miRNA [84], and induced cell cycle arrest [55,87]. miR-215 which had pro-apoptotic function was detected at high levels in normal human colon tissue but at low levels in many human colon cancer samples. Once again, differences in expression of this miRNA in B. tabaci Q and B might contribute to biological differences in developmental time, reproduction, and survival.

Conclusions
High-throughput sequencing enabled the study of miRNAs in B. tabaci, which is an important pest worldwide. We identified 170 conserved miRNAs and 15 novel miRNA candidates in B and Q. We found significant differences in the expression of miRNAs between B and Q, which might contribute to the displacement of B by Q. To date, little is known about the functions of these miRNAs in insects, especially in B. tabaci. Further analysis of the expression and function of these miRNAs could increase our understanding of regulatory networks in the insect and lead to novel approaches to its control.

Author Contributions
Conceived and designed the experiments: DC QG. Performed the experiments: QG YLT. Analyzed the data: QG. Contributed reagents/ materials/analysis tools: DC. Wrote the paper: QG DC.