Comparative transcriptome analysis of Gossypium hirsutum L. in response to sap sucking insects: aphid and whitefly

Background Cotton (Gossypium hirsutum L.) is a major fiber crop that is grown worldwide; it faces extensive damage from sap-sucking insects, including aphids and whiteflies. Genome-wide transcriptome analysis was performed to understand the molecular details of interaction between Gossypium hirsutum L. and sap-sucking pests, namely Aphis gossypii (Aphid) and Bemisia tabacci (Whiteflies). Roche’s GS-Titanium was used to sequence transcriptomes of cotton infested with aphids and whiteflies for 2 h and 24 h. Results A total of 100935 contigs were produced with an average length of 529 bp after an assembly in all five selected conditions. The Blastn of the non-redundant (nr) cotton EST database resulted in the identification of 580 novel contigs in the cotton plant. It should be noted that in spite of minimal physical damage caused by the sap-sucking insects, they can change the gene expression of plants in 2 h of infestation; further change in gene expression due to whiteflies is quicker than due to aphids. The impact of the whitefly 24 h after infestation was more or less similar to that of the aphid 2 h after infestation. Aphids and whiteflies affect many genes that are regulated by various phytohormones and in response to microbial infection, indicating the involvement of complex crosstalk between these pathways. The KOBAS analysis of differentially regulated transcripts in response to aphids and whiteflies indicated that both the insects induce the metabolism of amino acids biosynthesis specially in case of whiteflies infestation at later phase. Further we also observed that expression of transcript related to photosynthesis specially carbon fixation were significantly influenced by infestation of Aphids and Whiteflies. Conclusions A comparison of different transcriptomes leads to the identification of differentially and temporally regulated transcripts in response to infestation by aphids and whiteflies. Most of these differentially expressed contigs were related to genes involved in biotic, abiotic stresses and enzymatic activities related to hydrolases, transferases, and kinases. The expression of some marker genes such as the overexpressors of cationic peroxidase 3, lipoxygenase I, TGA2, and non-specific lipase, which are involved in phytohormonal-mediated plant resistance development, was suppressed after infestation by aphids and whiteflies, indicating that insects suppressed plant resistance in order to facilitate their infestation. We also concluded that cotton shares several pathways such as phagosomes, RNA transport, and amino acid metabolism with Arabidopsis in response to the infestation by aphids and whiteflies.


Background
Plants are sessile organisms; hence, they are easy target for biotic and abiotic stresses, as they cannot escape from these stresses. One million or more phytophagus insect species use plants as a source of their food. These phytophagus insects obtain nutrition from plants either by chewing or by sucking sap from aerial or underground plant parts. During their long period, approximately 350 million years of co-evolution with plants [1], both insects and plants evolved a variety of different interactions. These interactions can be either positive interactions (insect-mediated pollination and seed dispersion) or negative interactions (insects used plant parts as food). In contrast to the fitness of insects to use plants as a food source, plants have also evolved distinct mechanisms that deal with these interactions [2]. Specialized defense mechanisms in plants protect them from insects. The mechanisms are either inherent constitutive defenses, which include physical barriers such as cell wall and cuticle, or induced defense mechanisms in response to insect attacks [3]. The induced defense mechanisms also include the activation of proteinase inhibitor [4], polyphenol oxidase, chitinases, and so on [5] and the release of secondary metabolites, which attract the parasitoid of attacked insects [3]. The induced defense mechanism has recently been demonstrated to be mediated through SA (Salicylic acid) [6], JA (Jasmonic acid) [7], or ET (Ethylene) pathway [8]. Further, several reports support the fact that there is complex crosstalk between plant hormonal pathways that control the plant responses to wounds, insects, and pathogen attacks [9]. However, detailed genetic regulatory mechanisms that govern plants' interactions with insects are yet to be fully understood. Sap-sucking insects cause severe damage to both crop plants and glass house-grown plants, in both temperate and tropical regions [10]. The damage caused by the sap sucking pests can be direct as well as indirect. Direct damage occurs due to the removal of nutrients by the feeding of insects on plant phloem, which reduces plant vigor and causes shoot and leaf distortion. A more serious problem is due to the large amount of honeydew that they secrete onto leaves and fruits. The honeydew is colonized by sooty molds. The sooty mold interferes with photosynthesis and may lower the quality of fruits and vegetable harvest. Whiteflies and aphids serve as vectors for plant viruses such as Gemini virus, resulting in the spread of viral diseases [9]. The whitefly Bemisia tabaci is known to attack more than 500 species of plants, representing 74 plant families. They are particularly serious threats to crops such as squash, melons, cucumbers, pumpkins, tomatoes, eggplant, potatoes, cotton, and okra. Thus, recently, there has been a keen interest in studying the molecular interaction between sap-sucking insects and plants using microarray studies [3,9,10]. These studies identified the involvement of not only defense-related metabolism but also the genes related to normal cell metabolism such as cell wall modification, water transport, vitamin biosynthesis, photosynthesis, carbon assimilation, and nitrogen and carbon metabolism during aphid attack in different plants such as Arabidopsis thaliana, Apium graveolens, and sorghum [10][11][12]. Similarly, the expression of trypsin proteinase inhibitor, lipoxigenase, and xyloglucan-endotransglycoxylase genes was shown to be upregulated, and rubisco subunit and ubiquitin carrier protein were down-regulated [13] during aphid attack. The difference in the performance of aphids and whiteflies on Arabidopsis thaliana mutant PAD4 indicate that besides some commonality of their interaction with plants, aphids and whiteflies also have a different and unique way of interacting with plants [3]. Cotton is a fiber and an oilyielding crop that is grown all over the world. Four species of cotton are usually grown worldwide [14]; however, the contribution of Gossypium hirsutum L. to the total lint cotton production is maximum worldwide [15]. The productivity of cotton is severely affected by both biotic and abiotic stress worldwide [16]. About 1326 species of insects have been reported to attack cotton plants worldwide [17] and among these species, aphids and whiteflies are one of the major pests for cotton agriculture [18]. Thus, we decided to study the response of cotton at a molecular level in response to the infestation by whiteflies (Bemisia tabaci Gennadius) and aphids (Aphis gossypii G.) by using Roche's GS-Titanium pyrosequencing.

Development of common data set for digital expression analysis
The digital expression analyses of differentially expressed contigs were performed by generating a common data set. All the assembled contigs of different events were again assembled in a single pipeline by using Roche GSassembler. A total of 14810 contigs were generated with an average length of 684.43 bp (Table 2 and Additional  file 3). The transcript per million (TPM) value for each transcript was calculated and normalized against the expression of housekeeping genes (see Methods). The total assembled transcriptome represents 84.7% contigs from C, 82.6% contigs from A2, 82.9% contigs from A24, 86.6% contigs from W2, and 79.6% contigs from W24. The quantitative profiling of transcriptome using DEGseq (R-bioconductor) revealed that 158, 465, 123, and 100 contigs were up-regulated in A2, W2, A24, and W24, respectively; whereas 876, 753, 1013, and 1048 contigs were found to be down-regulated as compared with the control in A2, W2, A24, and W24, respectively (Additional file 4).
Transcriptional response of cotton toward the infestation of whiteflies was faster as compared with aphids The cotton plant responds to whiteflies by differentially expressing 1218 transcripts at 2 h of infestation, the number of differential transcripts decreases to 1148 at 24 h of whitefly infestation. In contrast to the number of differentially expressed transcripts in case of aphid infestation which were 1034 at 2 h and increased up to 1136 till 24 h of aphid infestation ( Figure 1). The amount of transcripts at W2 was statistically higher (p ≤ 0.002) than A2. Thus, the result indicates that the response of cotton plants toward whiteflies was rather fast, which gets decreased with time; whereas cotton plants respond slower to infestation by aphids and it was similar to response of whiteflies 24 h infestation. In order to find fold-change between aphid-and whiteflyinfluenced transcript, we checked the average inducibility of these differentially expressed genes and found that fold change was almost similar in all conditions of up regulation as well as 2 h of down regulated genes (Additional file 5); hence, the conclusion based on the number of differentially expressed genes would be meaningful. Further to confirm linearity in expression pattern among the three biological replicate, we have performed microarray with Affymetrix's Cotton Chip, the result indicates that the correlation between the replicates was as high as 0.97 to 0.98 in control (Additional file 6) and 0.87 to 0.94 in whitefly infested sample  (Additional file 7). Thus, the expression variation between replicates were very low and hence our transcriptome of single sample truly reflect the expression profiling of the event circumventing the need of biological replicates.

Major changes in transcriptomes in response to infestation by aphids and whiteflies
We analyzed top five genes that were either induced or repressed at 2 h and 24 h of infestation by either aphids or whiteflies ( Table 3). The top three genes induced at 2 h and 24 h were found to be common, namely, inositol oxygenase, phosphate translocator, and transketolase; further, phosphate translocator was found to be in the top-induced gene list of infestation by whiteflies at 2 h. Since phosphate translocator, including triose phosphate translocator, was reported to increase the source of carbon in the form of sugar [19], this finding seems to be an important way by which sap-sucking insects increase their sugar concentration in sap. The other two genes uniquely induced in A2 conditions were NADH dehydrogenase and asparginase. The higher expression of asparaginase also suggests the flow of nitrogen source into the sap. Divol et al. (2005) [11] reported that several genes which were involved in nitrate and sugar remobilization in celery, including glutamine synthase, were induced by aphid infestation. Our study also suggested that expression of the transcripts related to the cellular amino acid and nitrogenous metabolism were induced after aphid attack (Additional file 8). The other two highly expressed genes in A24 include RNA binding protein and bHLH type of transcription factor; the induction of these regulatory genes indicates transcriptional reprogramming induced by aphids. RNA binding gene was also induced during W2 infestation (Table 3).
In W2, genes such as protein kinase, β-xylosidase 1, and oxidoreductase (DMR 6) were highly induced. The β-xylosidase 1 (AtBxl1) was reported to be involved in secondary cell wall hemicellulose metabolism and plant development [20]. The initial stage of infestation by whiteflies seems to result in changes in cell wall modification by beta-xylosidase 1 (AtBxl1) and signal transduction through protein kinases cascade. During a later phase of infestation by whiteflies, that is, W24, genes such as threonine aldolase, ferric iron binding, 4hydroxyphenylpyruvate dioxygenase, hydrolase, and adenine phosphoribosyl transferase 1 were highly induced (Table 3). Further, the later stage of infestation by whiteflies also leads to changes in amino acid metabolism by changing the expression of threonin aldolase and 4hydroxyphenylpyruvate dioxygenase (Table 3). Similarly, hydrolases may also be involved in the hydrolysis of glucosinolate, which gives the cynate and nitril, and these products are reported to be toxic to the attacking herbivores [21]. The down-regulated gene list in A2 and W24 includes DNAJ heat shock N-terminal domain-containing protein. The gene amino acid transmembrane transporter and ethylene-forming gene 1-aminocyclopropane-1-carboxylate oxidase (ACO) were down-regulated in the A2 condition ( Table 3). Rerouting of amino acid transportation in response to aphid infestation was reported [22], and the suppression in the amino acid transmembrane transporter was probably linked to infestation by aphids in cotton. Genes such as chlorophyll-binding and copper ion-binding genes were also suppressed during the initial phase of infestation by aphids. Copper-binding genes were also suppressed in A24 (Table 3). Genes such as phosphoenolpyruvate carboxylase 3, DNA repair protein RAD23, and nodulin family protein were suppressed in A24. Phosphlipase D   alpha 2 gene, which is involved in wound response, ET, and ABA [23] signaling, was suppressed in the later phase of infestation by both insects. Genes such as serine-type endopeptidase and cellulose synthase 1 were suppressed in W2. Similarly, genes such as plasma membrane intrinsic protein 2A, myo-inositol-1-phostpate synthase 2, and translocon at the inner envelope membrane of chloroplasts 110 were down-regulated in W2 (Table 3). Genes such as plasma membrane intrinsic protein 2A and protease inhibitor were also downregulated in W24. Thus, the results clearly indicate that infestation by aphids and whiteflies influences changes in transcriptomes in cotton to promote their infestation; whereas the cotton responds to infestation by expressing certain genes or pathways to counteract the herbivorous behavior of these insects.

Expression of defense-related transcripts in response to infestation by aphids and whiteflies
Next, we examined the expression of genes reported to be involved in plant defense in response to various pathogens and insects (Additional file 9). We identified that several kinases were down-regulated in infestation by aphids and whiteflies, which includes Enhanced Disease Resistance 1 (EDR1), MAP kinase 6, MAP kinase 16, and cell wall-associated kinase 5. The roles of these kinases are very well reported in the literature in pathogen-induced plant immunity [24]. Interestingly, the involvement of MPK6-mediated phosphorylation and an increase in the stability of ACS (1-aminocyclopropane-1carboxylic acid synthase) leads to the production of a high level of pathogen elicitor-induced ET [25] response, which is reported. In our study, ACS 5 was suppressed at W24; MPK6 was suppressed in A2, A24, and W24 (Additional file 9); and ET-forming enzyme 1aminocyclopropane-1-carboxylate oxidase was suppressed in A2 (Table 3). Thus, the suppression of ASC5, MPK6, and 1-aminocyclopropane-1-carboxylate oxidase (ACO) during W24 and A2 indicates the insectmediated suppression of the ET pathway. Among the highly up-regulated genes are CBSX1 − cystathionine beta-synthase (CBS) family protein, which was reported to be involved in cysteine amino acid metabolism and also in cell redox maintenance [26]. Another gene that was highly up-regulated during these insect infestations was MLO1 representing mildew resistance locus 1; this gene was reported to be involved in the susceptibility of plants toward fungal pathogens [27]. Further, the involvement of callose synthase in response to pathogen attack and wounding was reported [28], and the involvement of glucan synthase in callose deposition and aphid resistance [29] was recently reported. We identified that a homologue of ATGSL10 representing glucan synthase was up-regulated in all the events; whereas homologues of ATGSL08 and ATGSL12 genes representing glucan synthase were down-regulated at most of the infestation time points, especially during the later phase of infestation (Additional file 9). The differential expression of callose synthase may be considered a strategy of compatible interaction that is used to prevent the plugging of the creation of pores by these insects [30]. We identified that senescence-associated gene 18 (SAG18) was upregulated in all the events (Additional file 9), indicating a probable strategy of channelizing the flow of free amino acids formed due to the breakdown of leaf proteins during senescence [31]. Further, the expression of Ca 2+ -binding DND2 [32] and CNGC2/DND1 was downregulated in all the events (Additional file 9); this may lead to an increase in the influx of free Ca 2+ to sieve elements to plugged damaged sieve elements due to sapsucking insects and prevent the loss of phloem sap [28].
Other genes representing the wall-associated kinase were highly induced during the initial phase of infestation. Similarly, in the differential expression of these cell walls homeostasis management related genes, other genes related to protein folding and cell membranes were also differentially expressed after infestation; for example, the expression of prefoldin 3 and 5 was up-regulated, whereas prefoldin 6 was down-regulated in infestation by insects. Further, membrane proteins such as SGR3, which represented Syntaxin/t-SNARE and SYP61syntaxin of 61 family proteins of plants, were downregulated in A2, A24, and W24 during infestation by these insects (Additional file 9). Interestingly, we identified that one of the major key proteins involved in microRNA biogenesis, that is, the AGO4-argonaute family protein, was found to be up-regulated in all the events, and the role of argonaute in plant immunity was recently reported [33]. Phenylalanine ammonia-lyase-2, the first enzyme involved in the phenylpropanoid pathway, was found to be up-regulated in W2 and W24. Similarly, the role of cytochrome P450 in the various biotic and abiotic stresses is well discussed [34]. We identified that the P450 genes family member CYP86A8 was induced in infestation by both aphids and whiteflies. The expression of the GTPase gene, which is also known as enhanced disease resistance 3, was induced in A24 and W2; whereas its expression was down-regulated in A2 and W24. Similarly, the expression of pathogenesisrelated 4 possessing chitinase binding activity was induced in A24 and W2; whereas its expression was suppressed in A2. Further, we identified the expression of NPR1-like protein 4, which was reported to be involved in plant resistance [35], was down-regulated in A2 and W24. Interestingly, a homologue of Arabidopsis NAC domain containing protein 2/ ATAF1, which is a negative regulator of plant resistance [36], was downregulated in the later phase of infestation by aphids and whiteflies. Plant WRKY DNA-binding transcription factors are involved in plant pathogen interactions [37]. We identified that the expression of WRKY 33 was enhanced in all the events (Additional file 9). However, the expression of WRKY 21 was down-regulated in A2 and W24; whereas WRKY 20 was down-regulated in A2, A24, and W24. WRKY 1 was down-regulated in A2, A24, and W2; WRKY 35 was down-regulated in W24; and WRKY 3 was up-regulated in all the cases. Thus, our results suggest that aphids and whiteflies interact with cotton plants via complex molecular interactions involving several pathogenesis-related genes and pathways.

Oxidative stress
Oxidative radicals play an important role in plants during various stresses, including the biotic stress such as insect infestation [38]. Thus, we checked the expression of genes that are involved in the scavenging of oxidative radicals (Table 4). Glutathione is a major ROS (Reactive oxygen species) scavenger in plants; we observed the expression of enzymes involved in glutathione synthesis, namely, glutathione peroxides in A2 and W24, ascorbate peroxidase (APX6) [39] in whitefly insect-infested leaves were down-regulated (Table 4). Further, catalases are the H 2 O 2 and other ROS detoxifying enzymes produced at the site of ROS/H 2 O 2 production [40]. We observed a decrease in the expression of catalase1 in all the four events, namely, A2, A24, W2, and W24; further, the expression of catalase 2 was decreased in W24, and a decrease in the level of catalases indicates an increase in H 2 O 2 level during the infestation by sap-sucking insects.
We further observed a decrease in the expression of superoxide dismutase family protein T5P19.1 in A24 and W24; copper chaperon for SOD (Superoxide dismutase) in all the events; copper/zinc superoxide dismutase3 (CSD3) in W24; copper/zinc superoxide  (Table 4). Thus, our results indicate that a decrease in the expression of the ROS scavenging enzyme may enable an increase in the concentration of ROS and H 2 O 2 , which are directly toxic to insects [38] and drive fast peroxidasemediated oxidative cross-linking of structural proteins in the cell wall [41].
GO classification and mining of the changes in the functional class after attacks by sap-sucking insects We used Arabidopsis thaliana model for 'GO' annotation of our transcriptome data. Differentially expressed genes of each event in comparison to the control were analyzed and functionally categorized based on three GO categories at p-values ≤ 0.05, which were performed by using singular enrichment analysis (SEA) of agriGO tool (http://bioinfo.cau.edu.cn/agriGO/). Only the molecular functions (F) and biological process (P)-related GO categories containing more than 4% genes of agriGO annotation were selected for further analysis (Additional file 8). These results showed that the differentially expressed transcripts were involved in various metabolic and cellular processes. The result showed that the major transcriptomic reprogramming in aphids happens in the late phase of infestation and most proportion of these categories were suppressed at later phase of infestation by both insects (Additional file 8). The major groups of up-regulated genes in aphid infestation belong to carboxylase, hydrolase, structural moleculer activity, and stress response by various signals and pathogens. Cellular catabolic and transporter related transcript were suppressed during aphid infestation. Further, the early up-regulated genes also shared the pathways belonging to a developmental or reproductive structure (Additional file 8). We observed that the later phase of infestation by aphids leads to the up-regulation of transcripts related to secondary metabolic processes such as phenyl propanoid biosynthesis, flavanoids, and aromatic compounds (Additional file 8). During the later phase of aphid infestation, leads to the suppression of transcript for phosphorus metabolic process, macromolecule metabolic process including post translational protein modification, RNA metabolic process, nitrogen compound metabolic process, transcript related to purin binding as well as ATP and adenyl nucleotide binding genes (Additional file 8). The result indicates the involvement of secondary metabolic pathways in infestation control of aphids, especially in the later phase. In addition, the later phase of aphid infestation (A24) also showed the up-regulation of genes belonging to amino acid and aromatic compound metabolism. We also noticed a decline in the expression of genes belonging to transporter activity in A2, and a significant number of genes belonging to pyrophosphatase and hydrolases of acid anhydrides are down-regulated in infestation by aphids; however, hydrolases of esters were up-regulated in A2. The transcriptomic reprogramming in response to whiteflies similar to infestation by aphids showed the up-regulation of several genes belonging to stress, response to signals and pathogens (Additional file 8). Interestingly, there was no major secondary metabolic pathway reprogramming in case of infestation by whiteflies in contrast to infestation by aphids (Additional file 8). Another contrasting difference was that many of the transporter activities were up-regulated in case of whiteflies (Additional file 8). Interestingly, a significant proportion of up-regulated genes belonged to transcription regulators, indicating that a later phase of infestation by whiteflies induces transcriptional reprogramming (Additional file 8).

The response of cotton plants toward aphids and whiteflies shares with hormonal and other biotic stresses
The members of genes responding to different hormonal pathways were obtained by querying transcriptome data to the Genevestigator database (https://www. genevestigator.com/gv/plant.jsp). We observed that transcripts responding to abscisic acid (ABA) were increased during the later phase of infestation by both aphids and whiteflies (Figure 2A). Transcripts responding to Jasmonic acid and salicylic acid were relatively higher in whitefly (W2) attack which was decreased with due course of infestation (W24). Similarly, transcripts responding to SA was higher during intial infestation of whiteflies (W2) which was in agreement with previous report where author suggests induction of SA pathways during the attack of whitefly on Arabidopsis [3]. We also observed that JA responsive transcript were suppressed during time course of whitefly infestation while there expression were enhanced in aphid infestation. These result support the previous report that inspite some commonalities in infestation mode of these insect, plant deals with them in different ways. JA-and SA-mediated induction of plant defense in response to insect infestation was indicated [3]. We have identified the expression of OPR3 (oxophytodienoate-reductase), which is involved in the JA biosynthesis, and the development of plant defense [42] was up-regulated in A24 and W2 but down-regulated in A2 and W24 (Additional file 9). The role of ACX (acyl-CoA oxidase) genes involved in JA biosynthesis was reported earlier [43]. We found that among the five family members of acyl-CoA oxidase (ACX), four of them were differentially expressed in our experiment (Additional file 9). ACX1 gene was upregulated in W2 and down-regulated in A24; ACX2 was up-regulated in A24, W2, and W24; ACX3 was downregulated in all except W2; and ACX4 was downregulated in all cases except W24. Similarly, expression of genes such as Jasmonate Resistant 1 (JAR1) was upregulated in W2 (Additional file 4). JAR1 encodes a JA amino acid synthetase that is involved in conjugating jasmonic acid to Ile. The result indicates the involvement of these pathways in the later phase of infestation management by plants. A decrease in the expression of lipoxygenase I (LOX1) and LOX2 in initial phase of infestation of these insects and a decrease in the expression of TGA2 in plants infested by aphids and whiteflies showed insect-mediated suppression of plant defense, which facilitated the compatible infestation of these insects (addional file 4). Transcripts related to ethylene and Giberlic acid was almost similar in all conditions but a significant increase is seen the the transcript related to W2 in cytokinin responses. We next examined the common genes with that of bacterial, viral or fungal infestation. Intrestingly response of cotton to both whitely and aphid show strinking similarity with bacterial response as compare to fungus and virus ( Figure 2B). The transcripts similar to bacterial response were increased in whiteflies during later phase of infestation (W24), similar trend was also observed in case of transcript similar to fungal response. The percentage of transcript similar to virual response were found to be almost similar in all conditions. The result showed differential expression of wound induced protein (WIN2), heat shock proteins (Additional file 4). Interestingly the transcript of ferritin is highly induced during whitefly infestation (Additional file 4). These insects also influenced the expression of RNAi machinarry. The ferritin is used to homeostatsis of Fe in cytoplasm of cell.
Recently, Kieu et al. [44] reported that Fe is necessory for proliferation of bacterial disease in Arabidopsis plant. The induction of ferritin after whiteflies infestation, raised interesting question wheather it is in favour of plant or insect, which need to be address in future. A decreased in the expression of non-specific lipase, which is involved in JA biosynthesis during infestation [45], shows the insect mediated suppression of JA pathway (Additional file 4). Similarly, an increase in the expression of downy mildew-resistant 6 (DMR6) (Additional file 10A), which is a positive regulator of susceptibility to fungi [46], also indicates a probable mechanism that induces a favorable response by insects in host plants (Additional file 4). Further, we have screened the pathogenic organism that also influenced the expression of these insect (aphid and whitefly) infestation-responsive transcripts by the help of Genvestigator (Additional file 11). Among them, fungi such as Alternaria brassicicola, Botrytis cinerea, Blumeria graminis, Erysiphe cichoracearum, E .orontii, Golovinomyces cichoracearum, Phytophthora infestans, and P. paraistica; bacteria such as Escherichia coli and Pseudomonas syringae; and viruses such as cabbage leaf curl virus (CalCUV) and turnip mosaic virus (TuMV) were found.

Different biological pathways operating during infestation by aphids and whiteflies
The differentially expressed transcripts compared with the control for each event were analyzed for the identification of biological pathways that were enriched using KOBAS (http://kobas.cbi.pku.edu.cn/home.do). During initial phase of infestation of these insects, citrate cycle (TCA cycle), lysine degradation, alpha-Linolenic acid metabolism, protein processing related transcripts were enriched in aphid infested plant while transcript of fructose, mannose, and sulfur metabolism were enriched in whitefly infested plants ( Table 5). Transcript of pathways like pentose phosphate, sesquiterpenoid and triterpenoid biosynthesis were enriched in both aphid and whitefly infested plant during initial phase (Table 5). Interestingly we obsereved that at later phase of whitefly infestation, transcript of amino acid metabolic pathways of glycine, serine, threonine, alanine, aspartate, glutamate, valine, leucine, isoleucine and histidine were enriched (Table 5).
We also observed that whitefly infestation also influence the expression pattern of transcript related to secondary metabolic pathway of flavonoid biosynthesis and vitamin metabolism specially ascorbate at later phase. Transcripts of fatty acid, porphyrin and chlorophyll metabolism were uniquely enriched in whitefly infested cotton plant while in aphid infested cotton plant transcript related to pantothenate and CoA biosynthesis, proteasome, galactose, arginine and proline, butanoate, endocytosis were enriched. Transcript of pathways related to peroxisome and phagosome were enriched in later phase of infestation of both insects. Transcript related to fatty acid elongation, circadian rhythm, brassinosteroid biosynthesis, branched dibasic acid metabolism were enriched during later phase of aphid infetstation. Expression pattern of transcript related to photosynthesis specially carbon fixation were significantly influenced by infestation of these insects ( Table 5). The transcripts of different pathways like carotenoid biosynthesis, glyoxylate and dicarboxylate metabolism were enriched in A2, W2 and W24 while beta-Alanine, pyruvate, and glutathione metabolism related transripts were enriched in A2, A24 and W24. Transcripts related to glycolysis and gluconeogenesis pathways were enriched at A2 and W24.
The differentially expressed transcriptomes of cotton in response to aphids and whiteflies share many commonalities with those of Arabidopsis thaliana L To further evaluate whether cotton and A. thaliana share a common pathway in response to infestation by aphids and whiteflies, we compared differentially expressed transcriptomes with those of publically available microarray data sets, namely, GEO: GSE6516 (whitefly induced) [3] and GEO: GSE5525 (aphid induced) [47]. We identified 798 and 134 common differentially expressed transcripts ( Figure 3A and 3B) between our transcriptome data and public domain microarray data for aphids and whiteflies, respectively. We further queried 798 and 134 genes against the KOBAS database to understand the common pathways shared between A. thaliana and cotton. Our result showed that expression of transcript related to carbohydrate metabolism, ribosome, RNA transport, phagosome, pyruvate, butanoate, and glyoxylate were commonly influnced in both arabidopsis and cotton plant during aphid infestation ( Figure 3C). In case of whiteflies, the number of common genes identified for A. thaliana and cotton were significantly less ( Figure 3B). The common genes in whiteflies exhibit the enrichment of glycerolipid metabolism, ascorbate and aldarate metabolism, glutathion metabolism, ubiquinone-terpenoid quinone biosynthesis, and protein processing in endoplasmic reticulum ( Figure 3D) in both A. thaliana and cotton. Aphids and whiteflies suck the sap from phloem; so, to evaluate whether there is any differential expression at the transcriptomics level in phloem cells, we compared the differentially expressed transcriptomes of aphid-and whitefly-infested data with publically available phloem transcripts (microarray data sets viz., GEO: GSE10247, Laser microdissected phloem cell-LMPC). We identified 190 and 212 common differentially expressed transcripts ( Figure 4A and 4B) between our transcriptome data and the public domain microarray data of the phloem-expressed transcripts for aphids and whiteflies, respectively. We further searched the pathways involved in the phloem cell after infestation by aphids and whiteflies and found that in whiteflies, sulfur metabolism and selenocompound metabolism-related transcripts were differentially expressed; whereas in aphids, oxidative phosphorylation-related transcripts were differentially expressed ( Figure 4C and 4D).

Validation of transcriptomic data using quantitative realtime PCR
Six differentially expressed genes including four upregulated and two down-regulated contigs from each condition as in the case of W2, W24, A2, and A24 were selected on the basis of their differential expression pattern in transcriptomic data for further validation. In case of W2, up-regulated genes like ADP-ribose pyrophosphohydrolase, hypothetical protein (contig00503), hypothetical protein (contig05537), pyrophosphatase and down-regulated genes like plasma membrane intrinsic protein 2A and inositol-3-phosphate synthase were selected. Similalry in case of W24, up-regulated genes like hypothetical protein (contig13398), NADPH dehydrogenase, oxidoreductase, trehalose-phosphatase and down-regulated genes like uridylyltransferase-related and hydrophobic protein (contig14199) were selected. In case of A2, up-regulated genes like hypothetical protein (contig05119), hypothetical protein (contig00504), cytoskeleton protein, hypothetical protein (contig15596) and down-regulated genes like ACC Oxidase and hypothetical protein (contigs 80) were selected. Similalry in case of A24, up- regulated genes like hydrophobic protein (contig02797), protein kinase, H 2 -translocating pyrophosphatase, ATPdependent peptidase and down-regulated genes like phospholipase D and DNA repair protein RAD23 were selected. The validation of these contigs for each condition, namely, C, W2, W24, A2, and A24, was carried out using three independent biological replicates. All the four induced contigs that were selected for W2 showed an expression that was 2 to 12 fold higher in the whitefly infested condition as compared with their respective control in qRT PCR ( Figure 5A). Similarly, in case of W24, the four induced contigs that were selected showed 2 to 65 fold higher expression as compared with non-infested controls ( Figure 5C). The down-regulated contigs of W2 and W24 showed their down expression in qRT PCR ( Figure 5B and D). Thus, the qRT PCR results on the contigs were selected in complete agreement with the transcriptional data. In case of aphids at both the events, namely, A2 and A24, the four induced contigs that were selected and represented both the events showed higher expression in the aphid-infested condition as compared with the non-infested control; however the range of expression varied from 1.5 to 2 fold ( Figure 5E and G). The down-regulated contigs of A2 and A24 showed their down expression in qRT PCR ( Figure 5F and H). Further, five contigs viz., 60S ribosomal protein L5, gene representing protein binding, kinase, 60S ribosomal protein L31 and EF-1-alpha were selected for real time validation, which showed constant expression in all the experiments, and this expression was observed as complementing the transcriptomes (Additional file 12). Thus, qRT PCR results agreed with transcriptomic data, however in case of A2 induced genes, the fold on induction over the control was relatively low ( Figure 5G).

Discussion
In the present study, we report the transcriptomic changes in Gossypium hirsutum L. leaf, in response to two sap-sucking insects (aphid and whitefly). Cotton plants were infested by these two insects, and transcriptome sequencing at an average 4.4X coverage was completed for the control and infested leaf samples. We observed that plants respond to whiteflies quickly by changing their transcriptome; whereas in case of aphids, the response is slow (Figure 1). The number of downregulated genes that were more than the up-regulated genes in infestation by both aphids and whiteflies support the previous report which showed that aphids stimulate the suppression of more genes than does induction [47]. Our study suggested that aphids and whiteflies influence the expression of cell and cell wall metabolism by changing the expression of enzymes of sugar metabolism such as phosphoenolpyruvate carboxylase 3, sugar translocator/posphate translocator, cell wall modifier β-Xylosidase 1, inositol oxygenase and cellulose synthase 1 (Table 3). We also identified that amino acid metabolism was significantly altered by changing the transcription of key enzymes such as threonin aldolase and 4-hydroxyphenylpyruvate dioxygenase (Table 3). These insects reroute the amino acid transportation [22], and cotton plants probably respond to it by the suppression of the amino acid trans membrane transporter, as a defense strategy that is deployed by plants in response to infestation by aphids. In case of sap-sucking insects, the amino acid composition of plant sap determines the attractiveness of insects [48], and sapsucking insects, thus, become the secondary sink of amino acid for plants and increase the genes related to the amino acid biosynthesis pathway [10]. It was also reported earlier in A. thaliana that amino acid permease AAP6 mutant reduces the amino acid level in phloem sap and this correlated with the aphid behavior [49]. The infestation mediates the up-regulation of senescence in response to aphids and indicates the breakdown of leaf proteins and probably the translocation of the free amino acid pool, thus forming the phloem sap [31]. The result of cotton transcriptome in response to infestation by aphids and whiteflies showed significantly enrichment of the amino acid biosynthesis pathway (Table 5). During infestation, these insects damage the sieve tube; in response to this damage, plant respond by the release of Ca 2+ , which causes plugging of the sieve plate [28] and prevents the loss of phloem sap. However, aphids overcome this defense by secreting Ca 2+ -binding protein through their saliva, thus preventing clogging. In our experiment, suppression of the Ca 2+ -binding protein was noticed, and this may be considered an insect-influenced plant strategy for increasing the Ca 2+ level in phloem sap by suppressing these genes (Additional file 9). In qualitative terms, the impact of 2-h infestations by aphids was similar to 24-h infestations by whiteflies; like cytokinin, fungus, bacteria, viruses, JA, gibberellin, and cytokinin responsive transcripts were similarly induced or suppressed by both conditions (Figure 2A and B). Further, we have screened the pathogenic organism that also influenced the expression of these insect (aphid and whitefly) infestation-responsive transcripts with the help of Genvestigator (Additional file 11). Among them,  fungi such as Alternaria brassicicola, Botrytis cinerea, Blumeria graminis, Erysiphe cichoracearum, E. orontii, Golovinomyces cichoracearum, Phytophthora infestans, and P. paraistica; bacteria such as Escherichia coli and Pseudomonas syringae; and viruses such as cabbage leaf curl virus (CalCUV) and turnip mosaic virus (TuMV) were found. In contrast to chewing insects, weak wounds were created by these phloem feeders. The artof-style insertion of these insects may be comparable to the fungus haustoria and bacterial infection response. It was reported earlier that intercellular fungal hyphae growth resembles with that of style penetration of whiteflies [50]. The GO annotation of the differentially expressed genes for A2, A24, W2, and W24 showed the involvement of various metabolic and cellular processes (Additional file 8) during infestation by these insects. The transcriptomic reprogramming in response to infestation by aphids and whiteflies showed the upregulation of several genes belonging to stress, response to signals and pathogens (Additional file 8). Some of the interesting features include the differential expression of transporters in response to sap-sucking insect infestation (Additional file 8); some of transporter-related transcripts are up-regulated in W2, whereas they are downregulated in A2. The inducibility of water transporter [11], gluatathion S conjugate transporter [22], and sugar transporter [50] was also reported earlier in response to infestation by aphids. In W2, the genes related to various developmental processes such as seed development, post-embryonic development, multicellular development, and root development were found to be induced (Additional file 13). The relationship between developmental genes and biotic stress was reported earlier; for example, seed development genes were found to be induced in response to nematode infestation [51]. The transcripts belonging to secondary metabolic processes such as phenyl propanoid biosynthesis, flavanoids, and aromatic compounds (Additional file 8) were upregulated during the later phase of infestation by aphids, and these aromatic compounds may be involved in the attraction of parasitoid of aphids [52]. The transcripts of some of the hydrolyses and carboxylesterase were enriched in 2 h of infestation with aphids (Additional file 8). The differential expression of hydrolase and transferase in response to biotic and abiotic stress was shown in the form of the differential expression of glycosyl hydrolase family 1 in P. rapae infestation in Brassica oleracea [53] and UDP-glycosyltransferase activity in toxic detoxification during insect infestation [54]. The hydrolysis of the product of glucosinolate gives the cyanides and nitriles, which are toxic to herbivores, and the transportation of these to the phloem leads to the deterrence of herbivorous and phloem-feeding insects. The role of glucosinolates in plant defense is well reported; they form in plant tissue and are transported to the site of insect attack [55]. In our result, enzyme benzoate-CoA ligase, which is involved in benzoyloxyglucosinolate synthesis, and genes such as nitrate transporter 1.9 and transporter protein containing the properties of transporters of glucosinolate [56] (Additional file 13) were differentially expressed in response to sap-sucking insect infestation. Similarly, the induction of glutathione Stransferases members in response to cabbage aphids [57], cell wall modification enzymes pectin transferase [58] was reported earlier. The sap-sucking insects are the chief mediators or vectors of spreading plant viruses.   There is a decrease in the expression of Poly (A) binding protein 2 (PABP2) in both insect-infested leaves, which help in the transmission of sap-sucking, insect-mediated virus infection and the translation of viral RNA molecule [59], heat shock protein 70 [60], and chloroplast RNAbinding protein 29, which are used in viral protein folding (Additional file 13). Salicylic acid (SA), jasmonic acid (JA), and ethylene signaling pathways are involved in the regulation of plant-induced defense after attacks by pathogens and insects [61]. The JA-responsive pathway is usually activated when there is an attack by necrotrophic and chewing insects; whereas the SA mediates the defense response against biotrophic pathogens and insects such as aphids and whiteflies [62]. Aphids and whiteflies influenced the expression of cotton JA and ET synthesis genes in our experiment. The role of JAR1 in JA-mediated defense development has been already reported. Further, the constitutive expression of JA-and ET-signaling pathways in CEV1 (Constitutive expression of VSP 1) mutant of Arabidopsis, which was resistant to aphid growth, was reported [62]. We identified that at a later phase, the infestation of whiteflies leads to the suppression of CEV1 expression (Additional file 4). Thus, our results indicate the fine tuning of the JA pathway in cotton in response to the infestation by aphids and whiteflies. Further, in addition to changes in the expression of the genes involved in JA and ET biosynthesis, we also identified that the expression of hormonal signaling kinases, including MAP2K9 and MAPK6, was also altered, and the relation of MAP kinases in defense mechanisms involving JA, SA, and ET is well established [63]. We also identified that the expression of enzymes involved in oxidative radical scavenging were suppressed after the infestation by aphid and whiteflies (Table 4); these may lead to an increase in oxidative radicals and H 2 O 2 in the phloem sap, which is a probable strategy that is deployed by cotton plants against insect infestation [38]. We also report the involvement of ABA and GA pathways during the infestation by aphids and whiteflies in cotton (Figure 2A). The role of ABA [64,65] and GA [66] in plant-insect interaction has been recently shown. Further, the involvement of ABA and GA during the defense responses against green bug phloem feeding in sorghum [67] has been recently demonstrated. It has been reported that the increase in GA causes the increase in trichome density, and this may protect the plants from aphids and whiteflies [66]. Our results further showed the involvement of cytokinin in defense responses to aphids and especially whiteflies (Figure 2A). The involvement of cytokinins in defense responses toward aphids or whiteflies has not been reported earlier; however, the role of cytokinins in plants and bacterial interactions [68] and bacterial isopentenyl transferase (ipt) genes, which are involved in cytokinin biosynthesis, which, in turn, are involved in resistance to the tobacco hornworm and green peach aphid nymphs, was reported [69]. The down-regulation of genes such as overexpressors of cationic peroxidase 3-OCP3 (Additional file 9 and Additional file 10B), non-specific lipase [45], LOX1 [70], and TGA2 [71] which leads to the suppression of phytohormonal-mediated plant resistance and increase in the expression of DMR6 (Additional file 4 and Additional file 10A), which is a positive regulator of the susceptibility of plants to pathogens [46], showed insectmediated suppression of plant defense and compatible infestation of these insects. Plants assimilate heavy metals such as Ni [72], Zn [73], and Se [74] for protection against herbivorous insects. We identified the enrichment of selenometallo metabolism in the case of whiteflies infested both cotton and Arabidopsis plant ( Figure 4D). Selenium is a member of sulfur(s) group, and, hence, plants readily assimilate selenate in place of sulfur into cysteine as selenocysteine (SeCys) via the sulfur metabolic pathway [74]; this explains the enrichment of the sellanometallo metabolic pathway in response to whiteflies. We also observed the enrichment of the transcript related to RNA transport both in cotton and Arabidopsis plants during aphid infestation ( Figure 3C). It is already reported that viruses hijack the plant RNA transportation system for disease spreading. In parellal, plant activate the pathway of mRNA surveillance to control formation of aberrant RNA, which is a defense mechanism, was also generated in response to virus infection in plants [75]. Sap-sucking insects are potential vectors for plant-borne viruses [10]; thus, our result also proposes an interesting question as to whether plants understand the potential threat of virus infection after the infestation by aphids and whiteflies.

Conclusions
An average of 4X coverage transcriptome information will be helpful in understanding the induced defense responses operating against aphids and whiteflies in agriculturally important cotton plants, and will also pave the way for developing new insect pest-management strategies. The expression pattern of transcripts reveals that sap-sucking insects interact with plants by suppressing the expression of positive regulators of phytohormonalinduced resistance genes, inducing the negative regulator of the plant resistance gene and suppressing the defense-related transcription factors such as WRKY and other MAP kinases involved in plant defense. Our result also suggests that these insects shift the sucrose and amino acid mobilization by changing the expression pattern of different genes related to amino acid and carbohydrate metabolism. Thus, insects facilitate their infestations and plants try to repel these insects by activating their glucosinolates, secondary metabolite pathway, and reactive oxygen weapons.

Plant material and insect infestation
The seeds of Gossypium hirsutum (var.MCU5) were sown for germination in a mixture of solarite, vermiculite, garden soil, and sand soil (1:1:1:1 ratio). After germination of seeds, plants were grown for five weeks in a glass house at 28 ± 2°C (day/night), a relative humidity of 50-60%, a 14-h photoperiod, and a photosynthetically active radiation of 900 μmol m-2 s-1 , respectively. Fiveweek-old plants containing six leaves were selected for the experiments. Cultures of aphid and whiteflies were maintained in potted cotton plants in the laboratory at 26 ± 2°C and 70% relative humidity [76]. About 25 freshly emerged whiteflies and 2nd instars aphid nymphs were released an average per leaf of plants. Whole experimental plants were covered with perforated polyethylene bags to prevent the insects from escaping. The insect infestation experiments were performed in three biological replicates. Solutions of 10 mM MgCl 2 (pH 7.0 /dissolve in PBS) were spread as mock solution [77]. Insects were removed by a fine brush after 2 and 24 h of infestation, and immediately, two middle leaves were frozen in liquid N 2 for total RNA isolation. All the experiments with Aphids and whiteflies were performed with approval of IBSC(Institutional Biosafety Committee).

RNA isolation and sample preparation
Total RNA were isolated by using Spectrum plant total RNA Kit (Sigma-Aldrich) according to the manufacturer's protocol and underwent DNaseI treatment. Out of three biological replicates, only one plant's RNA was used for transcriptome sequencing. To amplify the mRNA, double-stranded cDNA were prepared using oligo-dT primers containing T7 promoter and SuperScript W Double-Stranded cDNA Synthesis Kit (Invitrogen). These double-stranded cDNA were amplified using an in vitro amplification system of the Gene chip IVT labeling kit (Affymetrix). Amplified cRNA underwent double-stranded cDNA preparation by using random hexamer primer and SuperScript W Double-Stranded cDNA Synthesis Kit (Invitrogen). These double-stranded cDNA were purified by the QIAquick PCR purification column (Qiagen). The double-stranded cDNA (3 μg) were used for transcriptome sequencing. The transcriptome sequencing was performed as per the manufacturer's protocol of Roche's GS-Titanium pyrosequencing. Further to confirm linearity in expression pattern among the three biological replicate, microarray experiments were performed with isolated RNA sample through Affymetrix kit as per the manufacturer's protocol.

Assembly and annotation of transcriptomes
The reads from each library were assembled separately by Roche Newbler (GS-assembler) version 2.3. The assembly criteria were set as 40 bp overlap size and 90% identity between the reads. Contigs of each library were pooled together and assembled to form a common data set. The reads were counted from their respective libraries for the newly assembled contigs, and their TPM (Transcript per Million) was calculated. For further analysis, TPM values were log2 transformed (Additional file 14). Genes such as Actin (AY305733), UBQ7 (DQ116441), Gbpolyubiquitin-1 (AY375335), Gbpolyu biquitin-2 (EE592464), Histone 3 (AF024716), and 18S rRNA (L24145) [78] were selected for the normalization of the expressed contigs in each condition. Digital gene Expression was carried out by DEGseq package in Rbioconductor 2.15 for each library with reference to control. Contigs with p-value ≥ 0.05 were selected for differential gene expression. Two-fold up-and two-fold down-regulated contigs were selected. Therefore, 8 categories were made (control vs. aphid 2 h up, control vs. aphid 2 h down, control vs. aphid 24 h up, control vs. aphid 24 h down, control vs. whitefly 2 h up, control vs. whitefly 2 h down, control vs. whitefly 24 h up, and control vs. whitefly 24 h down) (Additional file 4). The differential genes for each transcripts were subjected to 2 by 2 Chi-square test. The test was performed between the differential genes of A2 vs W2 and A24 vs W24. The p-value for A2 and W2 was found to be 0.002 showing the 99% significance level while in case of A24 and W24 we got the p-value of 0.809 only. The contigs of each event were subjected to blast using program blastx with the TAIR 9 protein database (Additional file 15) and blastn for cotton ESTs available in the NCBI database (Additional file 16) at e-value 10 -5 .

Functional annotation
The TAIR IDs of the contigs in each event were used for the GO annotation. The detailed GO annotations were studied using the agriGO tool (http://bioinfo.cau.edu.cn/ agriGO/), which was categorized in biological processes, molecular functions. The differential genes were querid against the hormonal and biotic stress related transcripts in genevestigator (https://www.genevestigator.com/gv/ plant.jsp). All the differentially expressed genes were also subjected to KOBAS analysis (http://kobas.cbi.pku.edu. cn/home.do), and significant pathways were selected at the p value ≤ 0.05. Differentially expressed genes were also compared with the public databases generated from plants of Arabidopsis thaliana that were infested with aphids and whiteflies at different time points (GEO: GSE5525, GEO: GSE6516) and Laser Microdetection Phloem Cells (LMPC -GEO: GSE10247), which were derivatives of Arabidopsis thaliana. The genes that were