The Immune Strategy and Stress Response of the Mediterranean Species of the Bemisia tabaci Complex to an Orally Delivered Bacterial Pathogen

Background The whitefly, Bemisia tabaci, a notorious agricultural pest, has complex relationships with diverse microbes. The interactions of the whitefly with entomopathogens as well as its endosymbionts have received great attention, because of their potential importance in developing novel whitefly control technologies. To this end, a comprehensive understanding on the whitefly defense system is needed to further decipher those interactions. Methodology/Principal Findings We conducted a comprehensive investigation of the whitefly's defense responses to infection, via oral ingestion, of the pathogen, Pseudomonas aeruginosa, using RNA-seq technology. Compared to uninfected whiteflies, 6 and 24 hours post-infected whiteflies showed 1,348 and 1,888 differentially expressed genes, respectively. Functional analysis of the differentially expressed genes revealed that the mitogen associated protein kinase (MAPK) pathway was activated after P. aeruginosa infection. Three knottin-like antimicrobial peptide genes and several components of the humoral and cellular immune responses were also activated, indicating that key immune elements recognized in other insect species are also important for the response of B. tabaci to pathogens. Our data also suggest that intestinal stem cell mediated epithelium renewal might be an important component of the whitefly's defense against oral bacterial infection. In addition, we show stress responses to be an essential component of the defense system. Conclusions/Significance We identified for the first time the key immune-response elements utilized by B. tabaci against bacterial infection. This study provides a framework for future research into the complex interactions between whiteflies and microbes.


Introduction
Insects interact and coexist with various types of microorganisms in many different ways and have evolved sophisticated strategies both to recognize and degrade entomopathogens, as well as to benefit from bacterial mutualists [1,2]. Over 20% of all insect species are known to be associated with endosymbionts [3]. Amongst these, phloem-sap-feeding insects such as whiteflies, psyllids and aphids have been shown to possess specialized bacteriocytes that harbor primary and secondary endosymbionts [2]. These endosymbionts provide essential amino acids, protect the host from pathogen infection and also help it to adapt to different environments [3,4,5]. How these insects protect themselves from bacteria pathogens, while retaining beneficial endosymbionts, however, remains to be discovered [6]. Elucidating the host defense mechanisms of these insects will not only shed light on this question, but will also facilitate the development of novel insect-pest control strategies.
The whitefly, Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae), is a cryptic species complex composed of at least 36 morphologically indistinguishable species [7,8,9,10]. This species complex has members that rank as some of the most economically damaging insect pests [11]. Bemisia tabaci damages plants by direct sucking and by transmitting plant viruses [12]. In nature, whiteflies interact with various bacterial species, some of which are entomopathogenic and may serve as potential bio-control agents [13,14]. Previous pyrosequencing analysis has revealed a diverse range of bacteria present in B. tabaci, though latest report showed the diversity of bacterial communities in B. tabaci is relatively limited [15,16]. On the other hand, the primary (or obligatory) endosymbiont Portiera is considered to play a role in the synthesis of essential amino acids and carotenoids to the whitefly host [17,18], while secondary endosymbionts may regulate the life parameters in various ways [19,20]. These different interactions provide an excellent opportunity to study the immune system of an insect species that interacts simultaneously both with bacterial pathogens and endosymbionts. However, the molecular basis of the interactions between whitefly and those microbes remains largely unknown.
A previous study examined the transcriptome of whitefly exposed to the entomopathogenic fungus Beauveria bassiana and showed that only a limited number of canonical immune related genes were involved in the host's defense [21]. In addition, a genome-wide analysis and functional study showed that another hemipteran, the pea aphid, seems to lack many genes that are essential for the immune response in many other insects [22]. To investigate the immune system of Hemiptera further, we examined the gene expression profile of the whitefly B. tabaci under challenge from a well characterized Gram-negative bacterium Pseudomonas aeruginosa [23,24,25].
Most of the previous knowledge on insect immune response is based on host reactions after injection of bacteria into the insect's body cavity [1]. For most insects, however, the normal route of bacterial invasion is via oral ingestion [26,27]. More recent work has shifted from cavity injection to systematic intestinal immune responses under oral infection, which mimics the natural mode of bacterial invasion [28]. Here, we examined the whitefly's essential host-defense strategies after oral bacterial challenge. Functional analysis of the differentially expressed genes indicate that MAPK cascade, antimicrobial peptides (AMP) and gut epithelium renewal play critical roles in the whitefly defense system, whereas stress responses are also induced to improve host tolerance.

Plants, whiteflies and bacteria cultures
The Mediterranean (MED) species of the B. tabaci complex was used in all the experiments [7]. A culture of MED was maintained on cotton plants (Gossypium hirsutum L. cv. Zhemian 1793) in climate chambers at 2761uC, 14 h light/10 h darkness with 70610% relative humidity (RH). The purity of the MED colony was checked every 5 generations by RAPDs and with sequencing of the mitochondrial cytochrome oxidase 1 gene [29], [30]. The Pseudomonas aeruginosa strain ATCC9027 was obtained from the Microorganisms Germplasm Bank of Guangzhou, China.

Bioassay
Bioassays were carried out at 2761uC and 70610% RH. The infection solution was obtained from an overnight bacterial culture. The density of bacteria was adjusted to 1610 8 CFU/ml in 10% sucrose solution. The same sucrose solution without P. aeruginosa was used as control. Feeding and body cavity injection are two main approaches to carry out infection. In this study, the former method was chosen to mimic a natural mode of bacterial infection as well as to avoid physical damages to the whitefly. Transparent plastics tubes (L10, ø4 cm) were prepared as the feeding chambers for whiteflies. One end of the chamber was covered with a sandwich of 2 layers of carefully stretched Parafilm membrane separated by a layer of 1 ml bacterial solution. Approximately 100 newly emerged adult whiteflies were fed for 60 hours in each tube. The dead whiteflies were counted and cleaned out of the container every 6 hours. All treatments were replicated four times.

Sample preparation for sequencing
Approximately 1,500 newly emerged adult whiteflies were collected for each treatment. At first, the control and 24 h treatment groups were fed, respectively, with sucrose and bacterial solution. After 18 hours, the 6 h treatment group was fed with bacteria. At 24 hours post-infection (hpi), approximately 1,000 whiteflies were collected from the control, 6 and 24 hpi treatments, respectively. This method ensured that the whiteflies were collected at the same developmental stage. In addition, because some whiteflies may not feed on the artificial diet at the beginning, the 6 hpi and 24 hpi whiteflies were fed with bacterial solution throughout the treatment to make sure that every whitefly individual can ingest enough bacteria. Samples were frozen immediately in liquid nitrogen and homogenated using the FastPrep system (MP Biomedicals). Total RNA was purified with SV total RNA isolation kit (Promega) according to the manufacturer's instructions. RNA quality was assessed by Nanodrop 2000 (Thermo Scientific) and 2100 Bioanalyzer (Agilent) as previously described [21,31]. For each treatment, two biological replicates were conducted and processed independently. One replicate was used in the digital gene expression (DGE) library preparation and the other was used for real time quantitative PCR (qPCR) analysis.

Digital gene expression (DGE) sequencing and tag annotation
The methodology for DGE sequencing was largely based on that described in previous studies [31,32]. In brief, the mRNA from each sample was purified with magnetic oligo (dT) and subjected to cDNA synthesis. The cDNA was subsequently digested with NlaIII, which recognizes the CATG sites. Then adapter 1 was ligated to the site of NlaIII cleavage. The purified cDNA fragments were digested using MmeI that cuts 17 bp downstream of the CATG site, thus producing tags with adapter 1. Then, the adapter 2 was ligated at the site of MmeI cleavage. After 15 cycles of PCR linear amplification, 6% TBE polyacrylamide gel electrophoresis was used to purify the tags. After digestion, single strand molecules were added to the Illumina sequencing flowcell and fixed. The purified tags were sequenced by using Illumina HiSeq 2000 platform at the Beijing Genomics Institute (Shenzhen, China). DGE library data sets obtained from this work are available at the NCBI Gene Expression Omnibus under the accession number of GSE52837.
Clean tags were generated after removing 39 adaptor sequences, low quality sequences, empty reads and tags with a copy number of 1. A reference database containing all possible CATG+17 nucleotide tag sequences were created for the transcriptome of the MED whitefly (unpublished data, available upon request). Sequencing tags were mapped to the whitefly transcriptome reference database with no more than one nucleotide mismatch. The number of unambiguous clean tags for each gene was calculated for gene expression analysis and TPM (number of transcripts per million tags) was used to normalize the data.

Analysis of differentially expressed genes
The levels of gene expression were compared between: 1) the control library and the 6 hpi library; and 2) the control library and the 24 hpi library. False discovery rate (FDR),0.05 and the absolute value of log 2 Ratio$1 were used as the threshold to judge the significance of gene expression difference [33,34]. Gene Ontology (GO) classification system provides a dynamic, controlled vocabulary for all eukaryotes [35] and was used to annotate the possible functions of differentially expressed genes (DEGs). Also, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was used to depict the pattern of host response against bacterial challenge [36]. The number of DEGs in each GO term and KEGG pathway were also calculated. Using the MED transcriptome database as background, significantly enriched GO and KEGG pathway terms were determined using the hypergeometric test (P-value,0.05).

Real time quantitative PCR (qPCR) analysis
To confirm the results of the DGE analysis, the expression of 20 selected genes was measured using qPCR. cDNA was synthesized using the SYBR PrimeScript reverse transcription-PCR kit II (Takara). qPCRs were performed in 96-well plates using the ABI Prism 7500 fast real-time PCR system (Applied Biosystems) with SYBR green detection. Each gene was analyzed in triplicate, after which the average threshold cycle (C T ) was calculated per sample. The relative expression levels were calculated using the 2 2DDCt = 2 2[DCt (treatment)2 DCt (control)] method. The reference gene actin was used to normalize the expression level of other genes [37]. All of the designed primers were synthesized at Boshang BioCompany (Table S1).

Whitefly survival curve
This research focused on the response of whitefly to orally delivered bacteria; therefore it is important to find the important time points during bacterial infection. To achieve this, the survival rate of the whitefly was monitored for 60 h after P. aeruginosa ingestion. In this assay, a sharp decrease in the number of surviving whiteflies was observed during 12-36 hours postinfection (hpi) (Fig. 1). Therefore, 24 hpi was chosen as a sample collection time point. In addition, as sample collected earlier may reflect the whitefly's intestinal response to bacterial infection, the sample at 6 hpi was also collected for sequencing [26,38].

DGE library construction, sequencing and mapping
Three whitefly DGE libraries (Control, 6 hpi and 24 hpi) were sequenced and approximately 6 million raw tags were obtained for each library (Table S2). To determine the appropriate total tag sequencing number, a saturation analysis was carried out to check whether or not the number of detected genes kept increasing when the sequencing amount increased. Results for all three samples showed that when sequencing reached 5 million or higher, the increase in the number of detected genes was negligible (data not shown). After the removal of low quality reads, the numbers of distinct tags were 128641, 138701 and 126216 in the libraries of the control, 6 hpi and 24 hpi, respectively (Table S2). The ratio of clean tag to total tag was about 97% in each library. For annotation, the short tags of these DGE libraries were mapped to the MED whitefly transcriptome reference database. Tags which can map to more than one gene were filtered out. About 80% of clean tags were mapped unambiguously to the transcriptome database, showing the high quality of the sequencing and reference database. As a result, each library generated ,17,000 tag-mapped genes, which also means that about 36% of genes in the whitefly transcriptome could be detected during DGE sequencing (Table  S2).

Differentially expressed genes (DEGs) and qPCR validation
After annotation, we compared the control library with 6 hpi and 24 hpi bacterial challenged libraries to identify the DEGs that may play a central role in the host's defenses. Compared with the control, 948 genes were up-regulated and 400 genes were downregulated at 6 hpi, while 758 genes were up-regulated and 1030 genes were down-regulated at 24 hpi ( Fig. 2A, Table S3). Interestingly, at 6 hpi, the majority of DEGs were up-regulated. At this time point, the mortality of infected whiteflies started to increase, which indicated that the whiteflies had come into close interaction with the bacteria and had responded to pathogen infection. Therefore, 6 hpi reflects the early phase of whitefly defense against P. aeruginosa. At the 24 hpi, MED death increased rapidly. The bacteria had already clearly imposed substantial stress on the host. As a result, the sample at 6 hpi and 24 hpi represent different stages of MED's defense response.
The detected fold changes (log 2 ratio) of gene expression ranged from 29.57 to 9.57, and the majority of genes were up-or downregulated between 1.0-and 5.0-fold, respectively ( Figure 2B). In previous studies, the DGE method has been proven to have high reliability [31,39,40]. Due to the relatively high expense of Illumina sequencing, only one sequencing run was performed for each sample. To validate the DGE data, 20 selected genes were quantified for their transcription levels with qPCR in the 6 hpi and 24 hpi treatments. In 34 out of 40 tests, qRT-PCR results were consistent with the DGE data, providing further evidence of the reliability of our sequencing results (Table S1).  Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis are widely used to examine the biological processes in a large group of genes [35,36]. First, different GO and KEGG Orthology (KO) terms were assigned to DEGs. In the 6 hpi and 24 hpi data, 738 and 965 DEGs had GO annotations while 187 and 283 genes were annotated with KEGG terms. Enriched GO and KO terms in DEGs were then identified using hypergeometric analysis with the MED transcriptome database as background (Tables S4 and S5). In GO enrichment analysis, 124 and 136 GO terms under the category of biological process were over-presented in the 6 hpi and 24 hpi libraries,  (Table S5). More importantly, the phagocytosis and melanogenesis related pathways were enriched in DEGs suggesting the critical role of the whitefly's cellular immune responses during bacterial infection. The involvement of multiple GO and KEGG pathways indicated that P. aeruginosa infection activated a number of cellular and molecular responses, which is discussed below.

Microbial recognition and signal transduction pathways
Pathogen recognition is the initial event of pathway activation and systematic responses. Host defense responses are initiated when microbial molecules such as lipopolysaccharides and peptidoglycans are detected by pattern-recognition receptors (PRRs) [1]. Among those receptors, c-type lectins are sugar binding proteins specifically binding to polysaccharide chains on the pathogen's surface [41]. In infected whiteflies, two c-type lectins were up-regulated in both 6 hpi and 24 hpi treatment groups (Table 1). C-type lectins have been shown to enhance hemocyte encapsulation ability in cellular immunity and activate prophenoloxidase in humoral immunity [42]. Thus, MED c-type lectins are probably able to regulate cellular and humoral immunity upon bacterial infection.
Activation of signaling pathways follows pathogen recognition. The MAPK pathway, which comprises the ERK, JNK and p38 mediated kinase cascade, is a conserved insect host-defense repertoire [43]. Dysfunctions of JNK and p38 resulted in hypersensitivity toward infection and stress [44,45]. Here, several essential MAPK components were activated, including MAP3K7/ TAK1, its binding protein TAB1 and MAP3K13/LZK (Table 1). TAK1 is required for the activation of NF-kB and JNK pathways in Drosophila, while LZK is able to phosphorylate and activate JNK [46,47,48]. More interestingly, a putative peroxiredoxin was suppressed at 24 hpi and its homolog in Drosophila was identified as a negative regulator in the Tak1-JNK arm of the immune signaling system (Table S3) [49]. In addition, other genes associated with MAPK were activated, such as an insulin/growth factor receptor, a ribosomal protein S6 kinase and tribbles homolog 2. The findings provide additional evidence of the importance of the JNK pathways following P. aeruginosa infection. Moreover, the CREB binding protein/P300 and phosphoinositide 3-kinase, two players in JAK/STAT pathways, were also upregulated significantly ( Table 1).

Activation of AMPs and other effectors in immunity
Fast and massive production of AMPs has evolved in insects to be a central strategy of their immune system. Knottins are small proteins that have antimicrobial peptide activities and they are widely present in both plants and insects [50]. The DGE analysis revealed that three antimicrobial knottin genes (Btk 1, 2, and 3) were induced at 6 hpi and this was confirmed subsequently by the qPCR data (Table 2, Fig 3). In Drosophila, TAK1 and its binding protein are able to activate the NF-kB pathway and, ultimately, Table 3. Genes involved in cell proliferation and related pathway.  the production of AMPs [51,52]. Whether this signaling pathway is conserved in B. tabaci, however, remains to be discovered. Other than activation of AMPs, coagulation and melanization are also important for microbe sequestration and degradation. In insects, serine proteases and their inhibitors (serpin) are responsible for the activation and regulation of coagulation and melanization [53,54]. Several serine proteases are modulated in P. aeruginosa infected whiteflies. Interestingly though, 1 and 2 serine proteases are upand down-regulated at 6 hpi, whereas 2 and 6 serine proteases are clearly up-and down-regulated at 24 hpi, respectively ( Table 2).

Regulations of cell proliferation and epithelium renewal upon infection
The intestinal epithelium is the first protective barrier of the host from orally delivered microorganism infection. In Drosophila, stem cell division was activated upon infection to compensate for the cell damage due to the presence of bacteria and also to maintain gut homeostasis [38,55,56,57,58]. The enrichment of GO terms such as cell cycle (GO:0007049), positive regulation of cell proliferation (GO:0008284) and epithelium development (GO:0060429) indicated the involvement of cell proliferation and epithelium renewal in host defenses (Table S4). About 80% of genes were induced at both 6 hpi and 24 hpi and extensive genetic studies on the Drosophila gut have shown that this process is governed by JNK, JAK-STAT, Wingless and Epidermal growth factor receptor (Egfr) pathways [55,56,59,60], and all of these pathways were also regulated in our data. In addition, a Wingless like protein and several related genes in the Egfr pathway were induced (Table 3 were also activated, including SMAD1, SMAD4 and an insulin receptor (Table 3).

Stress response genes of whitefly
The high motility caused by P. aeruginosa clearly showed that a hyper-biotic stress had been imposed on the whitefly by 24 hpi. DEGs were consequently highly enriched in the GO term: regulation of response to stress (GO: 0080134) at both 6 hpi and 24 hpi (Table S4). This result is consistent with the activation of the MAPK, because their roles in the stress response are widely accepted [43].
Chaperones work in protein folding and quality control, which maintains cellular homeostasis and buffers environmental stress. Our analysis showed that several chaperones and cofactors were regulated. Five out of seven genes were induced at 6 hpi whereas five out of six genes were repressed at 24 hpi (Table 4). Induced genes include several cofactors of HSP90 and HSP40 proteins. In addition, a large group of genes involved in detoxification, such as Cytochrome P450, Glutathione S-transferase and Glutathione peroxidase were also regulated ( Table 5). These genes are likely involved in the detoxification of reactive oxygen species, as well as other toxins produced by bacteria [61,62]. Interestingly, 15 out 16 Cytochrome P450 genes were suppressed at 24 hpi. Cytochrome P450 genes were also suppressed following pathogen infection [56,63]. It suggests there is a general trend of P450 suppression after pathogen infection in other insect species.
Another notable event in infected whitefly was the activation of DNA repair proteins. Five and seven DNA repair genes were induced in 6 hpi and 24 hpi whitefly, respectively ( Table 4). The only repressed gene was a negative regulator in DNA damage repair (Table 4) [64]. DNA damage caused by pathogenesis might have triggered the DNA repair response, because several genes in mismatch repair and double-strand break response have been reported to be highly activated to improve host defense [65]. An alternative explanation, however, is that these genes might also participate in cell-cycle checkpoint, as our data also show that cell proliferation is activated [66].

Modulation of basal metabolism
Orally-delivered infection of P. aeruginosa also modulated a large group of basal metabolism processes, especially in the late infection stage. Given that intestinal function is disturbed in pathogenesis, these regulations may act as a self-protective strategy to maintain and reallocate energy supply. Our data shows that the majority of metabolism-related genes were down-regulated at 24 hpi (Fig 4). Suspension of feeding has been discovered in other insects that encountered bacterial infection and thus the inhibition of metabolism might be a direct result of this [67]. Whitefly may have evolved this adaptive strategy to prevent further pathogen ingestion. In addition, this down-regulation might be due to the disruption of energy supply under bacterial infection.

Discussion
Despite its economic importance as a global pest, little is known about the immune response of the B. tabaci complex to bacterial infection. Recent developments in next-generation RNA-seq technology, however, allow a systematic study of the whitefly's host defense strategies upon intestinal bacterial infection. In whiteflies collected at both time points following bacterial infection, more than one thousand genes were modulated.  Functional analysis uncovered the complexity of the host's responses to pathogenic bacteria. Bacterial infection not only induced genes in immune signaling and several types of effectors, but also altered the expression pattern of several sets of genes in xenobiotics detoxification, protein folding, DNA repair and basal metabolism. Altogether, our analysis provides the first rough outline of the whitefly's defense mechanisms employed against pathogenic bacteria and raises further questions that deserve more investigation.
The major goal of our research was to decipher the immune strategies of B. tabaci. Humoral and cellular immunity are the two arms of the insect's defense system, and in the former response antimicrobial peptides (AMPs), enzymatic cascades and other soluble effectors are employed to degrade foreign invaders [1]. In our study, both DGE data and qPCR analysis showed the activation of AMP production, highlighting again its core status in the antimicrobial response. As with the cellular immune system, pathway analysis showed clearly the involvement of phagocytosis and melanogenesis, more research is needed to construct a comprehensive image. In addition, our data implicates the function of intestinal stem cell proliferation in gut homeostasis maintenance, which may be conserved in whitefly.
The signaling pathways governing immune processes were also revealed by our work. We observed significant activation of several arms of MAPK cascades. Raf and ribosomal protein S6 kinase are essential kinases in the ERK pathway, while TAK1, TAB1 and TRB2 can control both p38 and JNK pathways [43,68]. Studies in Drosophila showed TAK1 and TAB are located at the cross point of the JNK and Imd pathway in fruit fly [51,52]. Upon activation, TAK1 and TAB can regulate NF-kB and then activate AMP production. Ras/MAPK signaling can also mediate intestinal homeostasis and regeneration in Drosophila. Exploring downstream events controlled by MAPKs, therefore, may provide key answers to how those processes are regulated in whiteflies. Besides the MAPK pathway, only a few genes in other canonical pathways were identified in infected whiteflies. For instance, we failed to identify the canonical factors in the Imd pathway from whiteflies. Interestingly, the pea aphid, another hemipteran insect, also appears to lack the Imd pathway [22]. Nevertheless, the absence of these genes may be due to the limitation of our reference database, which only accounts for a part of the B. tabaci genome. Likewise, although several antimicrobial knottins were induced, their induction level was relatively modest compared to that of other insects. Surprisingly, several defensins were not up-regulated upon pathogen infection (data not shown). More analysis is needed to characterize the function of whitefly AMPs.
Our analyses also enable a close examination of the stress response strategy of B. tabaci. Though stress response genes are not involved directly in immunity, their importance in the host's defense systems is recognized. Activation of these genes can help the host maintain cellular homeostasis and increase its capacity to endure the infection [69]. The involvement of chaperones, detoxification enzymes and DNA damage repair is likely to help the whitefly build up a high tolerance towards infection. In fact, these genes also participate in environmental adaptation and the development of resistance to insecticides [70,71]. An understanding of the regulation of these gene sets in particular may help in the development of novel insecticides.
In summary, we report for the first time the results of an NGS investigation into the molecular interactions induced by the oral delivery of a bacterial pathogen, P. aeruginosa, to the whitefly B. tabaci. Functional analyses of DEGs indicated that at 6 hpi both humoral and cellular responses are involved in the whitefly's defense responses. Furthermore, MAPK cascade, AMP and gut epithelium renewal probably play critical roles in the defense system, whereas stress response genes are also induced to build stronger host tolerance. Of particular interest is that only a few genes in other canonical pathways were identified in infected whiteflies. Further research built upon these findings may present an opportunity for the development of a novel whitefly control technologies.