Deep Sequencing-Based Transcriptome Analysis Reveals the Regulatory Mechanism of Bemisia tabaci (Hemiptera: Aleyrodidae) Nymph Parasitized by Encarsia sophia (Hymenoptera: Aphelinidae)

The whitefly Bemisia tabaci is a genetically diverse complex with multiple cryptic species, and some are the most destructive invasive pests of many ornamentals and crops worldwide. Encarsia sophia is an autoparasitoid wasp that demonstrated high efficiency as bio-control agent of whiteflies. However, the immune mechanism of B. tabaci parasitization by E. sophia is unknown. In order to investigate immune response of B. tabaci to E. Sophia parasitization, the transcriptome of E. sophia parasitized B. tabaci nymph was sequenced by Illumina sequencing. De novo assembly generated 393,063 unigenes with average length of 616 bp, in which 46,406 unigenes (15.8% of all unigenes) were successfully mapped. Parasitization by E. sophia had significant effects on the transcriptome profile of B. tabaci nymph. A total of 1482 genes were significantly differentially expressed, of which 852 genes were up-regulated and 630 genes were down-regulated. These genes were mainly involved in immune response, development, metabolism and host signaling pathways. At least 52 genes were found to be involved in the host immune response, 33 genes were involved in the development process, and 29 genes were involved in host metabolism. Taken together, the assembled and annotated transcriptome sequences provided a valuable genomic resource for further understanding the molecular mechanism of immune response of B. tabaci parasitization by E. sophia.


Introduction
The whitefly Bemisia tabaci (Hemiptera: Aleyrodidae), is well known as a worldwide invasive pest and may cause severe damage to various vegetables by feeding on phloem sap and transmitting many viruses [1]. It is a complex species containing at least 30 cryptic species [2]. B and Q-types are two most economically damaging and invasive species [3]. There are many studies focus on biological characterization, resistance, invasive mechanism, and biological control of B. tabaci [4][5][6][7][8][9][10][11][12]. Over the past years, B. tabaci has demonstrated a remarkable resistance to many groups of chemical insecticides [13][14][15][16]. Due to the rapid resistance development, it is necessary to explore an alternative and effective management strategy to control B. tabaci. Parasitoid or parasitoid-produced regulatory molecules can be used to improve conventional pest control strategies.
Endoparasitoids have been identified as very important natural enemies of various arthropods, and could be used as biological control agents [17][18][19]. Hymenopteran endoparasitoids deposit their eggs into the host insect haemocoel, whose larvae feed on the host until its death [20][21]. Encarsia sophia is one of the specific parasitoids of Aleyrodidae species and has been used as efficacious classic biological control agents in many regions [22]. It can parasitize all instar nymphs of B. tabaci, especially the third and fourth instar nymphs [23]. The female wasp is generated by a bisexual process, but the male wasp is produced by autoparasitism [24]. Homogeneous E. sophia prefers to lay male eggs in the host parasitized by the heterogeneous wasp. When E. sophia and other kinds of wasps are raised or released together, the antecedent colonizers should inhibit the colonization of followers [25]. Previous studies have shown that E. sophia has strong plasticity adaption abilities [26].
However, the relationships between endoparasitoids and their hosts are complicated and involve long-term co-evolution. Many studies have investigated parasitoid biological characteristics, chemical communication, phylogenetic co-evolution, and physiological responses [27]. An increasing number of researchers have focused on revealing the physiological mechanism underlying the parasite induced immune defensive system and the biological development of hosts in order to estimate the co-evolution process between parasitoids and their hosts [28][29][30][31]. Although several reports have concentrated on the molecular regulation mechanisms, there have only been a few descriptions of related, functional genes [32,33]. Furthermore, the limitations of previous research methods has led to the development of high-throughput RNA sequencing technology (RNA-Seq) [34].
RNA-Seq is widely used to obtain transcriptomes of the organism, tissue, or organ, to identify genes that were regulated under certain conditions, and to reveal the regulatory mechanisms in different organisms [35][36][37][38][39]. In recent years, RNA-Seq has increasingly being applied in the biological agents to reveal the interaction mechanisms in the complex parasitoid-host system. Transcriptome profiling of organism under parasitization helps us to obtain a better understanding of host responses and effect on host's growth, development. As a model species, Drosophila melanogaster and its parasitoid wasp Asobara tabida (Hymenoptera: Braconidae) is a well-studied system. Most genes associated with insect immunity appeared to be differentially expressed after wasp parasitized [40]. Most transcriptome studies on parasitoid-host systems have focused on Lepidoptera and Coleoptera, such as Plutella xylostella, Chilo suppressalis, Tenebrio molitor and Octodonta nipae [41][42][43][44]. A previous study showed that another parasitoid, Eretmocerus mundus may parasitize B. tabaci and induce the specific transcription of functional genes related to immune responses in the host [45]. However, the host manipulation by the parasitoid is species-specific, and the molecular mechanism of immune system in B. tabaci parasitization by E. sophia has not yet been explored. In this study, we used deep sequencing to explore B. tabaci response to E. sophia parasitization. Our results demonstrate that immuneand metabolic-related genes that are differentially expressed in parasitized versus non-parasitized B. tabaci nymph.

Insects Rearing and Parasitization
The biotype Q of Bemesia tabaci was obtained from the greenhouse at the Beijing Academy of Agriculture and Forestry. All experimental populations were derived from one pairs of newly emerged B. tabaci female and male. In our laboratory, the B. tabaci was reared on cotton plants (Zhong-mian-suo 49) in insect proof cages at 26 ± 1°C, and with a photoperiod of 15L: 9D. The purity of the cultures was monitored every three to five generations using the random amplified polymorphic DNA-polymerase chain reaction technique with COI gene [46]. E. sophia was obtained from the greenhouse at Beijing Academy of Agriculture and Forestry. All whitefly instar nymph stages were provided as hosts to E. sophia. Then approximate fifty E. sophia (female to male ratio of 8:1) individuals were released into cages to breed and the newly emerged female and male as parents for five generations breeding.
Thirty pairs of whiteflies were fed on cotton leaf in a micro insect cage and the fresh cotton leaf were provided every 24 hours. When they had reached later 3 rd or early 4 th instar, they were transferred in culture dish with a piece of cotton leaf, whose petiol were wrapped into soggy cotton, and then the mated E. sophia was released into B. tabaci rearing cage for parasitization. Sixty paired E. sophia were released into one culture dish. Wasp E. sophia were removed after 2 hours parasitization. The first group of samples was collected at 24-hr after parasitization (24AP). At this time period, the parasitoids were at the egg stage in which the embryo had formed and gradually began to move. The brown substance in the egg began to accumulate and chorion had appeared. In other words, the parasitoid possessed immune regulation ability, but the ability was not strong at the egg stage. Therefore, we could identify the immune defense response of the host against the parasitoid. The second sampling period was 72-hr after parasitization (72AP) when the wasps reach larval stage move around and absorb nutrition from the host. At this time, E. sophia may start to regulate host development and metabolism to finish their own development in whiteflies. Each treatment and control had three replicates.

cDNA Library Construction and Illumina Sequencing
Total RNA was extracted from all nymph samples using TRIzol™ reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instruction and treated with DNaseI. The concentration and integrity of RNA sample were determined using 2100 Bioanalyzer (Agilent Technologies). The first-and second-strand cDNA synthesis, end reparation, addition of "A" bases to 3' ends, ligation of adapters at the end of DNA fragments, and PCR amplification. The cDNA library was qualified and quantified with an Agilent 2100 Bioanalyzer and ABI StepOne-Plus Real-time PCR system, respectively, and then sequenced using the Illumina HiSeq™2000 platform at the Beijing Genomics Institute (BGI, Shenzhen, China).

Transcriptome Analysis
In order to obtain clean reads, the low quality and adapter-polluted reads were removed from raw data. The good quality reads were assembled using Trinity [47] and assembled sequences were output as unigenes. All raw sequencing data have been deposited in NCBI Sequence Read Archive (SRA) database (http://www.ncbi.nlm.nih.gov/sra) with the following accession numbers: SRR1909644 (24AP), SRR1909651(72AP), SRR1909652 (CK-24AP), and SRR1909653 (CK-72AP). All the open reading frames (ORF) of unigene in B. tabaci were identified. If a unigene had many ORFs, we selected the longest one.
The unigenes were used for BlastX search and annotation against the NCBI non-redundant (nr) (http://blast.ncbi.nlm.nih.gov/Blast.cgi), Swiss-Prot (http://expasy.org/tools/blast), Kyoto Encyclopedia of Genes and Genome (KEGG, http://www.genome.jp/kegg/) databases with an E-value cut-off of 10 -5 . Gene Ontology (GO) annotation of unigenes was analyzed using the Blast2Go software [48], and GO functional classification for all unigens was performed using the WEGO software [49]. In the absence of B. tabaci and E. sophina genome sequences, we selected eight transcriptome datasets of B. tabaci from the NCBI database, and try to utilize the annotation that were the most closely related to B. tabaci gene in the parasitized library.

Differentially Expressed Gene (DEG) Analysis
In order to find all the differentially expressed genes, the same FPKM (Fragments Per Kilobase per Million fragments) value of unigene was first calculated for the treatment and control groups [50]. The results were displayed as fold changes, p-values and q-values. According to the q-value (p-value's statistical result after PFR (Positive False Rate) correction), a q-value less than 0.05 or the absolute value of fold change greater than 2 represented a significant difference between the treatment and the control.

Quantitative Real-time PCR (qRT-PCR) Validation
The quantitative real-time PCR technique was used to verify the reliability of the deep sequencing. Nine differentially expressed genes were randomly selected. The β-actin gene was used for normalization. The four RNA samples represented nymphs at 24AP and 72AP, and their respective control (non-parasitized nymphs) at the same developmental stages.
First-strand cDNA was synthesized from the total RNA (1.2 μg) by using PrimeScripTM 1 st Strand cDNA Synthesis Kit (TaKaRa) with oligo (dT) 18 as primer following the manufacture's protocols. The reaction system consisted of 10 μl of SYBR Green, 0.4 μl of ROX, 2 μl of diluted cDNA, 0.4 μl of each primer and 6.8 μl of distilled water. The reactions were loaded on the CFX96™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA) under the following conditions: 50°C for 2 min; 95°C for 2 min; and 40 cycles of 95°C for 10s, 60°C for 15s, and 72°C for 20s, followed by melting curve generation (68°C to 95°C). Data analysis was performed by oneway ANOVA following by Tukey's test using SPSS software.

Illumina Sequencing and de novo Assembly
In order to know how E. sophia parasitization regulated B. tabaci development, immuneresponse and the differences in regulatory mechanisms between E. sophia egg-and larvaestages. Approximately, 35 million and 39 million reads were generated from non-parasitized and parasitized B. tabaci nymphs at 24AP, respectively, and 31million and 37 million reads were from non-parasitized and parasitized B. tabaci nymphs at 72AP, respectively. De novo assembly produced 292,696 B. tabaci unigenes with an average size of 616 bp. Of these unigenes, 35.96% were between 200 and 300bp, 27.43% were between 300 and 500bp, 22.65% were between 500 and 1000bp and 13.96% had nucleotide lengths above 1000bp (Fig 1).

Functional Annotation and Classification
For functional annotation, the 292,696 unigenes were aligned to the GenBank protein databases with a cut-off E-value of 10 −5 using BLASTx. Using this approach, 46,406 unigenes (15.8% of all unigenes) were successfully mapped. In order to predicate protein function, the unigenes were further given a gene ontology (GO) classification and subjected to KEGG pathway analysis. A total of 35,688 unigenes were annotated and assigned to GO terms, which consisted of three main categories: biological process, cellular component and molecular function.
A total of 11,993 unigenes were categorized as cellular components, 12,102 unigenes were grouped under the molecular function, and 11,593 unigenes under biological processes. KEGG pathway analysis indicated that there were 4,721 unigenes assigned to different pathways in which translation, signal transduction, neurodegenerative diseases, infectious diseases, and endocrine system were the main B. tabaci pathways after E. sophia parasitizzation.

Enrichment Analysis of DEGs
A total of 1,482 genes appeared to be significantly differentially expressed in the parasitized and non-parasitized B. tabaci, of which 852 genes were differentially up-regulated and 630 genes were differentially down-regulated (Fig 2A). At 24AP, there were 584 genes differentially expressed, of which 356 genes were up-regulated and 228 genes were down-regulated. At 72AP, there were 1,270 genes differentially expressed, of which 698 genes were up-regulated and 572 genes were down-regulated. Out of all of regulated genes, 202 up-and 170 down-regulated genes were found at both time points ( Fig 2B) and more genes were up-regulated than that of the down-regulated genes at both 24AP and 72AP (Fig 2A). Furthermore, there was a significant difference in the numbers of differentially expressed genes at 24 hours than at 72 hours after parasitization. When E. sophia emerge in the larvae stage, more genes seemed to be involved in regulatory responses as compared to the egg stage. During the larvae stage, the parasitoid could move freely and began to feed on the host tissues. The distribution of the regulated genes indicated that their expression levels (>95%) were between two-to six-fold higher than at the egg stage (24 AP). Only a few genes changed more than six-fold (Fig 2).
GO analysis revealed that the DEGs were mainly categorized in the cellular component cluster, that focus on macromolecular, organelle, and cellular levels. In the molecular function cluster, the DEGs were mainly found in structural molecule, binding, and catalytic activity. In the biological process cluster, the DEGs were mainly categorized in cellular and metabolic processes, and cellular component organization or biogenesis (Fig 3). In addition, more genes were involved in cellular processes, metabolic processes, single-organism processes, response to stimuli, biological regulation, localization, and cellular component organization or biogenesis at 72AP. Translation and signal transduction were the two most important pathways according to the KEGG pathways analysis. For KEGG enrichment analysis, genes involved in the immune system, nervous system, endocrine system, and metabolic activities were differentially expressed. The above results showed that parasitization had a great impact on the normal life activities of the host.

Effects of Parasitism on the Transcription of Host Immune-related Genes
Vertebrates have a set of immune defense mechanisms that include innate immunity and adaptive immunity, but invertebrates only have innate immunity protection [51]. Insects will initiate their innate immune response when encounting foreign agents, such as bacteria, fungi, virus, and parasitoid. The immune system of insects can be divided into two categories: 1) humoral defense, including the antimicrobial peptides, reactive intermediates of oxygen, melanin formation and clotting; and 2) cellular defense mainly based on haemocytes, such as phagocytosis, encapsulation, microaggregation and nodulation [52][53][54]. Two defense mechanisms are associated with a wide range of immune-related genes.
Our sequencing results indicated that E. sophia parasitism had a significant impact on the transcription of immune-related genes in B. Tabaci nymph (Table 1). We identified several upregulated genes with homologs known to be involved in immune responses in insects, such as: defensin, knottin, serpin I2, laminin, spectrin, and apolipophorin. Defensin is an antimicrobial peptide, which acts as an innate immunity effector molecule and provides the first protection from pathogen infection. After parasitization by E. sophia at 24AP, we found that the transcription of defensin was up-regulated in B. tabaci nymph. Our results were consistent with previous studies that the mRNA levels of defensin in D. melanogaster and Phlebotomus duboscqi were significantly increased after parasitization [55,56]. Although the main action targets of defensin are bacteria and fungi, it also plays a role in the host-parasitoid system. Knottins are mini proteins that are present in many different organisms and have various biological functions [57]. After parasitization by E. sophia at both 24AP and 72AP, four knottins were overexpressed. Like defensin, it is also an important antimicrobial peptide.
Serpin I2 was one of the genes having higher levels of up-regulation (5.33-fold) at 72AP in RNA-seq analysis. Quantitative RT-PCR analysis (Fig 4) also show that it was up-regulated by 6.85 folds. Serine proteases are important immune regulatory proteins which play a significant role in the activation of the prophenoloxidase (PPO) cascade. The cascade activation eventually causes melanization to kill parasitized wasp through choking [58], however, serine protease inhibitor (serpin) can prevent the serine proteases activated melanization and weaken host defense for wasp parasization. Although studies have shown that serpins can be regulated by the parasitoids infestation in many hosts, their transcriptional levels are different in different parasitoid-host systems, and even in the same parasitoid-host system, two opposite situations may occur. Mahadav et al. and Song et al. found that serpins were down-  regulated in parasitized B. tabaci nymphs and P. xylostella larvae [26,59], while Etebari et al. [41] discovered that serpins were up-regulated 2-to 7-fold after P. xylostella parasitization by Diadegma semiclausum. In C. chilonis parasitized C. suppressalis, three up-regulated and three down-regulated serpins were identified in the fatbody [42]. Different serpins may play different roles in immune defense. Cellular immunity is another important component of the insect immune system. Laminin can stimulate cell adhesion and cell movement. Cofilin is an actin-binding protein which promotes cell migration and movement by changing the adhesion between cells and the extracellular matrix. Actin plays a significant role in facilitating cellular activities. The up-regulation of these genes showed that the host enhanced hemocyte encapsulation by reinforcing the extension and adhesion of hemocytes. laminin, cofilin, and actin were identified in our study and they were over-expressed at 24AP and 72AP. Ras3 and Rho1 are related to cellular immunity in D. melanogaster [40,60]. In our study, two genes were also identified significantly differentially expressed after parasization which may also be involved in the immune response of B. tabaci. At 24AP and 72AP, these genes were consistently over-expressed, which indicated that the cellular immunity not only defend parasitoid embryo and larval attacking. Superoxide dismutase (SOD), peroxidase (POD) and catalase (CAT) are three common enzymes in organisms. Organisms produce reactive oxygen species (ROS) under environmental stresses, which is cytotoxic to cells. However, the organism utilizes these protective enzymes to eliminate redundant ROS and protect themselves from damage [61]. When Trichoplusia ni is infected by baculoviruses, the expression of manganese superoxide dismutase (MnSOD) significantly reduces oxidative damage [62]. Zhu et al. also discovered that the transcriptional levels of Tenebrio molitor superoxide dismutases were up-regulated following bacterial infection or parasitization by Scleroderma guani [63]. In our study, the expression of SODC, PERO, and CATA were suppressed two-to four-fold at 72AP, but did not show significant changes at 24AP. At 72AP, the E. sophia reaches the larval stage and the damage to host becomes worse than that of egg stage. The decrease in the transcription of protective enzymes showed the parasitoid immune suppressive strategy.
In our study, some genes involved in insecticide resistance or detoxification were found to be differentially expressed under parasitization. Although these genes have no direct connection with defense against parasitoids attacking, they can be regarded as a stress response caused by parasitoid secretions. Takeda et al. [64] confirmed that the activities of glutathione-S-transferase (GST) and cytochrome P450 (CYP) increased in parasitized P. xylostella larvae. We also found most of the cytochrome P450 genes were highly expressed after being parasitized by E. sophia. Some genes were over-expressed at 24AP and others were over-expressed at 72AP.
Heat shock proteins (HSPs) are recognized as a family of highly conserved chaperones which respond to all kinds of environmental stress factors, such as heat, toxins, UV radiation, and invading pathogens by protecting protein from misfolding and denaturation [65]. We identified three heat shock protein genes, which were homologous with D. melanogaster and Anopheles albimanus, and are involved in B. tabaci development. In addition, heat shock 70 kDa protein cognate 3 was found to participate in the immune response. Therefore, we deduced that heat shock protein families can defend the host from damage by participating in the immune response and B. tabaci development.

Effects of parasitism on the transcription of host development-related genes
The parasitoids complete their development by absorbing the host's hemolymph and tissues. However, the development of the parasitoid and the host are synchronous. A previous study found that Aphidiu servi parasitized Acyrthosiphon pisum late-stage nymph stopped growth [30]. B. tabaci late nymphs parasitized by Encarsia bimaculata also stop growing [66]. After Encarsia formosa parasitized Trialeurodes vaporariorum Westwood nymph, the wasp didn't molt to until host nymph reached to last instar [67]. In order to complete development, the parasitoids have to change the host's development to match their own growth. In some cases, parasitoids suppress host's development and accelerate the host's early-maturity [68], while, other parasitoids prolong host's development to meet their own developmental needs. Previous studies have proposed that the wasp might control host's development through regulating the juvenile hormone and ecdysone levels [69,70].
Juvenile hormone epoxide hydrolases (JHEHs) have been identified as regulatory proteins in the catabolism of juvenile hormones [71,72]. A previous study showed that JHEH transcript levels were down-regulated more than two-fold in P. xylostella after parasitization by D. semiclausum [41]. However, Wu et al. [42] discovered that JHEH and juvenile hormone esterase (JHE) transcript levels increased in C. suppressalis after C. chiilonis parasitization. Based on our transcriptome data, parasitization by E. sophia led to JHEH1 up-regulation at 72AP (Table 2 & Fig 4). However, up-regulation of larvae cuticle protein and down-regulation of pupal cuticle protein might imply that the parasitoid suppressed the host's development. Thus, the high concentrations of JH may lead to up-regulation of JHEHs and their activity in order to maintain the balance.

Effects of parasitism on the transcription of host metabolism-related genes
Stearoyl-CoA desaturase (SCD) is an endoplasmic reticulum enzyme that catalyzes the biosynthesis of monounsaturated FA from saturated FA [73]. SCD inactivation causes obesity and abnormal lipid metabolism and one SCD activity, SCD1, was induced by insulin, but inhibited by leptin [74]. We found that at 24AP and 72AP, the transcript levels of SCD in B. tabaci nymph were up-regulated 6.69 and 4.52 times, respectively (Table 3). Furthermore, genes involved in the insulin signaling pathway were also significantly up-regulated. Our result implied that the wasp regulated the lipid metabolism of the host in order to get more nutrients available in host and meet their own needs. A report showed that the wasp preferred to parasitize late instar larvae because of adequate nutrition [75]. Stearoyl-CoA desaturase is an essential enzyme for the parasitic Trypanosoma brucei, and RNA interference of SCD caused a reduction of the parasitemia and an increase in host survival [76]. Environmental stress can influence the organism's metabolism, same as parasitoid infestation, which is energetically consumption process [77]. We found a high number of differentially expressed transcripts were related to organism metabolism. Metabolic changes occurred at both time points, but a greater amount and different kinds of genes were affected at 72APthan at24AP. E.sophia infection influenced carbohydrate, lipid, and energy metabolism in the host. Some studies have found that trehalose content changed after parasitization [78,79]. In two treatment groups, maltose was degraded to glucose under the action of maltase. Beside the upregulation of maltase, other genes, in the citrate cycle and glycolysis, were over-expressed at 72AP, such as citrate synthase, aconitate hydratase, and glyceraldehydes 3-phosphate dehydrogenase. Glycolysis and the citrate cycle are carbohydrate metabolisms to produce ATP. Citrate synthase, aconitate hydratase and succinyl-CoA synthetase are three essential enzymes in the TCA cycle. Up-regulation of genes that control the synthesis of these enzymes showed that total ATP decreased in the organism. Therefore, the insect needed to obtain more energy by increasing the reaction rate of the TCA cycle. In addition, we found that the transcriptional level of cytochrome c oxidase and f-type H + -transporting ATPase were significantly enhanced. Cytochrome c oxidase is involved in ATP synthesis as a terminase of the mitochondrial inner membrane respiratory chain [80]. However, whether it can be regarded as evidence of enhanced respiration is not clear. Measurement of respiration rate should be investigated in future studies.
There are three types of ion transporting ATPases: P-type, V-type, and F-type. In organisms, their main function is to synthesize ATP and transport H + [81]. There were more overexpressed f-type H + -ATPases at 72AP than at 24AP. This suggests that E. sophia parasitization of B. tabaci involved increased energy consumption. The host was regulated to produce more energy to supply to the parasitoid. Visser et al. found most wasps lacked a lipid synthesis mechanism and could not accumulate energy [82]. Therefore, it is reasonable to assume that the parasitoid may continually obtain energy from the host in order to complete its development.

Conclusions
In summary, our study first presented comprehensive transcriptome profiles of B. tabaci in response to E. sophia parasitization using RNAseq. The most of differentially expressed genes of B. tabaci after parasization have potential roles in immunity, development and metabolism to meet parasitoids needs. The transcriptome profiles provided a basis for future research in elucidate the host-parasitoid interaction. In addition, the identified immune-, development and detoxification-related genes may be target for B. tabaci control.