Skip to main content

Transcriptome analysis of Gossypium hirsutum flower buds infested by cotton boll weevil (Anthonomus grandis) larvae



Cotton is a major fibre crop grown worldwide that suffers extensive damage from chewing insects, including the cotton boll weevil larvae (Anthonomus grandis). Transcriptome analysis was performed to understand the molecular interactions between Gossypium hirsutum L. and cotton boll weevil larvae. The Illumina HiSeq 2000 platform was used to sequence the transcriptome of cotton flower buds infested with boll weevil larvae.


The analysis generated a total of 327,489,418 sequence reads that were aligned to the G. hirsutum reference transcriptome. The total number of expressed genes was over 21,697 per sample with an average length of 1,063 bp. The DEGseq analysis identified 443 differentially expressed genes (DEG) in cotton flower buds infected with boll weevil larvae. Among them, 402 (90.7%) were up-regulated, 41 (9.3%) were down-regulated and 432 (97.5%) were identified as orthologues of A. thaliana genes using Blastx. Mapman analysis of DEG indicated that many genes were involved in the biotic stress response spanning a range of functions, from a gene encoding a receptor-like kinase to genes involved in triggering defensive responses such as MAPK, transcription factors (WRKY and ERF) and signalling by ethylene (ET) and jasmonic acid (JA) hormones. Furthermore, the spatial expression pattern of 32 of the genes responsive to boll weevil larvae feeding was determined by “in situ” qPCR analysis from RNA isolated from two flower structures, the stamen and the carpel, by laser microdissection (LMD).


A large number of cotton transcripts were significantly altered upon infestation by larvae. Among the changes in gene expression, we highlighted the transcription of receptors/sensors that recognise chitin or insect oral secretions; the altered regulation of transcripts encoding enzymes related to kinase cascades, transcription factors, Ca2+ influxes, and reactive oxygen species; and the modulation of transcripts encoding enzymes from phytohormone signalling pathways. These data will aid in the selection of target genes to genetically engineer cotton to control the cotton boll weevil.


Cotton is a fibre and oil-yielding crop grown worldwide. Four species of cotton are generally cultivated [1]; however, Gossypium hirsutum L. contributes the most to the total lint cotton production worldwide [2]. Cotton productivity is severely affected by both biotic and abiotic stresses [3]. Approximately 1326 species of insects have been reported to attack cotton plants. Among these species, the aphid (Aphis gossypii G.), the fall armyworm (Spodoptera frugiperda), the budworm (Heliothis virescens), the cotton bollworm (Helicoverpa armigera) and the cotton boll weevil (Anthonomus grandis) are the major pests affecting cotton culture [4]. The cotton boll weevil is undoubtedly the most devastating pest [5, 6]. The adult female feeds, ovoposits, and develops primarily in cotton flower buds and fruits. After hatching, the larvae remain within the reproductive structures and use them as a food sources and a protected habitat to complete their life cycle. The endophytic behaviour of the larvae has made these insects difficult to control with conventional insecticides and other cultural practices [7].

Plants have evolved elaborate defence systems to respond in a rapid and effective way to herbivorous insects. Numerous studies have revealed that, in addition to the constitutive defences comprised of trichomes, thick secondary cell walls, and toxic compounds, plants are equipped with inducible defences that can be grouped into indirect responses, such as the production of volatile odour blends to attract natural enemies of the attacking insect, and direct responses, such as the production of anti-digestive proteins, toxic secondary compounds, and enzymes that affect insect growth and development [8]. These two powerful defence systems evolved by plants during the long arms race with herbivores have enabled plants to survive [9].

An appropriate defence response to a biotic threat requires initial recognition. Herbivores or pathogens are recognised when conserved patterns of molecules, called herbivore- or pathogen-associated molecular patterns, HAMP or PAMP (such as chitin or flagellin) are detected by pattern recognition receptors (PRR) on the surface of the host plant cell, leading to HAMP-triggered immunity (HTI). Damage-associated molecular patterns (DAMP), which are endogenous molecules produced by the plant after infection, are also recognised by PRR to trigger defensive reactions [10, 11]. Following the recognition of an attacker, plants use different signalling cascades to reprogram their phenotype. These include an interconnected network of signal transduction pathways depending mainly on the small regulators jasmonic acid (JA), salicylic acid (SA), and ethylene (ET), concomitantly with Ca+2 ion fluxes, mitogen-activated protein kinases (MAPK), transcription factors and reactive oxygen species (ROS) [12].

Crop losses due to insects constitute one of the most significant constraints to the increase of global productivity and food production, estimated at 10-20% for major crops [13]. A better understanding of the diversity of plant responses to insect attacks and, in particular, the induced defences and their regulation, has generated interest in the scientific community to examine alternative strategies to protect plants and crops from insects pests by exploiting the endogenous resistance mechanisms exhibited by plants to most herbivorous insects. Recent transcriptomics studies of plants exposed to herbivory identified a central role for transcripts that can lead to the development of insect-resistant crops [4, 1418]. High-throughput sequencing of RNA using next-generation sequencing platforms (RNA-Seq) offers a variety of new possibilities such as the transcriptional profiling of organisms lacking sequence information [19], as well as the identification of novel loci, alternative splicing events [20], and sequence variation [21]. Thus, we decided to study the molecular responses of G. hirsutum flower buds to infestation by cotton boll weevil (A. grandis) larvae using the Illumina HiSeq™ 2000 platform. Furthermore, we combined this approach with laser microdissection (LMD) to isolate RNA from two different regions of tissue damaged by feeding larvae to evaluate the spatial expression pattern of the differentially expressed genes in response to feeding by cotton boll weevil larvae.


Analyses of RNA-seq data

To explore the response of G. hirsutum to tissue-chewing pests such as cotton boll weevil larvae, the flower bud transcriptome of infested plants was compared with control plants (non-inoculated flower buds). Two biological replicates for each condition were selected for transcriptome analyses with high-throughput parallel sequencing using HiSeq™ 2000, Illumina. We generated 327,489,418 sequence reads, and each sample was represented by at least 74 million reads, a tag density sufficient for quantitative analysis of gene expression (Table 1) [22, 23]. The correlation between the two biological replicates was high (0.99 for infested flower buds and 0.98 for control samples; data not shown).

Table 1 Summary of sequencing data output, statistical analysis of the reads obtained and mapping of the reads onto the cotton ( Gossypium hirsutum ) transcriptome

The sequence reads were aligned to the G. hirsutum reference transcriptome (cotton EST database) using the BWA package, an efficient engine when searching for perfect matches [24]. Among the total number of reads, 39–43.2% were confined to exons, and 67,538,707 were perfect matches (OMM) to the reference sequence (Table 1). The total number of expressed genes (contigs) was higher than 21,697 per sample (Table 1), and 21,561 expressed genes were common to all samples analysed (data not shown). The average length of contigs generated was 1,063 bp (Additional file 1).

Overview of the changes in gene expression in response to feeding by cotton boll weevil larvae

The quantitative profiling of the transcriptome using DEGseq (R-bioconductor) analysis identified 443 differentially expressed genes (DEG) in cotton flower buds infested with cotton boll weevil larvae: 402 of them (90.7%) were up-regulated, and 41 (9.3%) were down-regulated compared to the control (adjusted p-value ≤ 0.05, |log FC| ≥ 2.0). Among these DEG, 432 (97.5%) were identified as orthologues of A. thaliana genes by Blastx with an e-value of 10−5 (Additional file 2). To examine the range of genes involved in the response of cotton flower buds to inoculation with A. grandis larvae, the Blast2GO program was also used to confirm the annotation of differentially expressed transcripts (Additional file 2). Blast2GO software returned functions for 87.3% of the differentially expressed genes from species with greater Blast hit distributions that included Vitis vinifera, Glycine max, Populus trichocarpa, and Threobroma cacao. Of the DEG, 76% had the same functional annotation between Blastx with Arabidopsis and the Blast2GO analysis (Additional file 2).

To determine which genes and pathways were relevant responses to cotton boll weevil larvae feeding, a gene set enrichment analysis (GSEA) approach was used. Based on GSEA, a hypergeometric test (p-values ≤ 0.005) was applied to identify which cellular components (CC), molecular functions (MF) and biological processes (BP) were overrepresented in our list of DEG. GSEA revealed many DEG putatively associated with the plasma membrane and cell wall in the cellular component category (Additional file 3). Most of the genes encoding plasma membrane proteins are receptor-like kinases (RLK), which are known to be involved in the perception of pathogen-derived elicitors. All RLK genes are up-regulated in infested plants in comparison to control plants (Additional file 4). Among the genes associated with the cell wall, there are up-regulated transcripts encoding a disease resistance protein (LRR) and cell wall-modifying enzymes, such as pectin methylesterase (PME41) and endotransglycosylase/hydrolase proteins (XTH31, XTH32, XTH23 and TCH4) (Additional file 5). Many of the down-regulated genes encode cytosolic heat shock proteins such as HSP90 and HSP70 (Additional file 5). These molecular chaperones assist in folding newly synthesised proteins and also in several other biological and cellular processes, such as cell growth, development and signal transduction during abiotic and biotic stress [25].

With regard to molecular functions, the DEGs were mainly associated with sequence-specific DNA binding transcription factor activity (Table 2), calmodulin binding and calcium ion binding (Additional file 3).

Table 2 The expression of 88 Gossypium transcription factor (TF) genes in response to cotton boll weevil larvae feeding in flower buds

Furthermore, GSEA revealed significant biological processes altered in plants upon infection with cotton boll weevil larvae. We found that biological processes related to hormone biosynthesis and signalling, the response to organic substances, the regulation of biological processes, the detection of biotic stimuli, systemic acquired resistance, the respiratory burst involved in the defence response and innate immune responses were overrepresented by GSEA (Figure 1). Moreover, biological processes associated with defence against insects, such as the response to chitin, the regulation of the plant-type hypersensitive response (HR), the regulation of programmed cell death and death were also represented (Figure 1). One hundred thirty-four transcripts were found in the “response to chitin” category (GO:0010200), and 36 of them were also annotated with the gene ontology (GO) term “death” (GO:0016265) (Table 3). Of these, many transcripts encoded proteins involved in signal transduction such as receptor-like kinases, mitogen-activated protein kinases, calcium ion binding proteins and calmodulin-like proteins. In addition to genes associated with signal transduction, genes associated with hormone biosynthesis and the 26S proteasome pathway were annotated with the GO term “response to chitin and death” (Additional file 6).The DEGs were mapped using MapMan to generate a representative overview of the pathways affected (Figure 2). This analysis indicated the involvement of several genes in the biotic stress response including receptors recognising microbe-, pathogen-, herbivore- and damage-associated molecular patterns (MAMP, PAMP, HAMP and DAMP) as well as genes involved in triggering defensive responses such as transcription factors and the ET and JA signalling pathways (Figure 2).

Figure 1
figure 1

Distribution of differentially expressed genes (DEGs) (x-axis) into Gene Ontology (GO) categories (biological process) (y-axis) according to Gene Set Enrichment Analysis (GSEA). Only biological processes (BPs) discussed in the results are presented here. A complete list of BPs can be found in Additional file 6.

Table 3 Subset of differentially expressed genes (DEG) in cotton flower buds in response to feeding by cotton boll weevil larvae
Figure 2
figure 2

Representative overview of DEG involved in biotic stress response in Gossypium hirsutum 48 h after infection with cotton boll weevil larvae. Log2-fold change of gene expression (cotton boll weevil compared to mock-inoculated control) was analysed by MapMan software. Yellow squares represent up-regulated genes and blue squares represent down-regulated genes. The colour saturation indicates log fold change > 4 and < 4. The figure shows that the molecular recognition of pathogens and herbivores by plants to trigger a defence response requires initial recognition including the following: 1. Microbe-, pathogen- and damage-associated molecular patterns (MAMP, PAMP and DAMP) are recognised by pattern recognition receptors (PRR) and lead to PAMP-triggered immunity (PTI). 2. Oviposition-associated compounds are recognised by unknown receptors and trigger defensive responses. 3. Putative herbivore-associated molecular patterns (HAMP) are recognised by receptors and lead to herbivore-triggered immunity (HTI). 4. Wounding by herbivores leads to the release of DAMP and to wound-induced resistance (WIR).

Main processes or pathways affected by the response to infection with cotton boll weevil larvae

Several subsets of genes were identified, including protein kinases, Ca2+-binding sensor proteins (CaM), genes responsive to oxidative stress, transcription factors, cell-wall modification genes and phytohormone-responsive genes.

Genes associated with signal transduction: kinases, Ca2+-binding proteins and protein degradation

An appropriate defence response to a biotic threat initially requires recognition of the threat. Herbivore-associated molecular patterns such as chitin are detected by pattern recognition receptors (PRR) on the surface of the host plant cell, leading to HAMP-triggered immunity. Well-characterised plant PRR comprises leucine-rich repeat (LRR) receptor-like kinases (RLK). Members of this family of proteins have a common structure formed by an extracellular ligand-binding domain and a cytoplasmic kinase domain [26]. Eleven RLK/LRR are differentially expressed (up-regulated) and annotated with the GO term “response to chitin” and/or “death” (Table 3), including receptor-like cytoplasmic kinase (RIPK), receptor-like kinase localised on the plasma membrane (RLK1), serine/threonine-protein kinase-like domain (CCR4 and CCR3) and putative leucine rich repeat transmembrane protein kinase (EVR). Recognition of conserved PAMP by PRR triggers intracellular signalling via mitogen-activated protein kinase (MAPK, MAPKK, MAPKKK) cascades followed by the activation of WRKY transcription factors [27]. In our analysis, cotton boll weevil larvae affected the expression pattern of several mitogen-activated protein kinase genes. For instance, two mitogen activated protein kinases (MAPK) that encode MAPK3 in A. thaliana were transcriptionally up-regulated (contig2100 and contig2102). Similarly, a member of the MAP Kinase Kinase family (MKK9, contig14909) and a protein kinase kinase kinase (MAPKKK14, contig7892) that play a role in leaf senescence during pathogenesis by Alternaria blight in A. thaliana were also transcriptionally up-regulated [28]. The activation of MAPK through wounding is dependent on the direct binding of calmodulins (CaM) in a Ca2+ dependent manner. Mechanical wounding and insect attack trigger a transient increase of cytosolic Ca2+ levels within minutes. This Ca2+ oscillation acts as a mediator for stimulus responses to several Ca2+-binding proteins, such as calmodulins (CaM), calcineurin B-like proteins, and calcium-dependent protein kinases (CDPK) [29]. In our analysis, nine transcripts coding for calmodulin-like proteins, including EDA39, CML38, CML24 and ACA2, were present in the DEG related to chitin response. In Arabidopsis, these genes are involved in plant innate immunity signalling and function as sensors of the signal generated by Ca2+ influx into the cytosol (Table 3) [30].

Targeted protein degradation via the ubiquitin/26S-proteasome pathway is another important regulatory process. Among the thirteen G. hirsutum genes annotated to be involved in this process, four were also related to the chitin and/or death GO (Table 3).

Cotton boll weevil-induced transcription factor (TF) genes

Transcription factors serve as important regulators of biotic and abiotic stress responses by turning on or off the immune system during plant defence. Our results suggest that the gene expression reprogramming in response to the boll weevil larvae feeding depends on many transcription factors (Table 2). As many as 88 TFs, which are 20% of the DEGs, were differentially expressed after cotton boll weevil infection, with 87 up-regulated and only one down-regulated (HSF, Heat-Shock transcription factor). Among the TFs up-regulated by cotton boll weevil feeding, there are 23 ethylene response factors (AP2/ERE), 21 WRKY, 12 C2H2 zinc finger, nine helix–loop–helix (bHLH), six C3H zinc finger, six MYB, six NAC and five GRAS (Table 2). Members of the AP2/ERE and WRKY families were also overrepresented in the study reported by Libault et al. [31], in which the expression patterns of Arabidopsis transcription factors were analysed in response to purified chitooctaose. Similarly, Wei et al. [14] observed that among the differentially expressed genes in Arabidopsis leaves in response to diamond back moth feeding, there were 18 ethylene response factors (ERFs) and 10 WRKY. The ERF subfamily belongs to the APETALA2 (AP2)/ethylene-responsive-element-binding protein (EREBP) family, an FT family exclusive to plants. WRKY transcription factors in plants are one of the largest families of zinc finger transcription factors and modulate development as well as responses to abiotic stresses, wounding, pathogens, and herbivore attack [32]. WRKY TFs can act as both positive and negative regulators of plant defence pathways. The mechanisms activating WRKY TF might involve the MAP kinase cascades and/or calcium signalling [3234].

In the present study, 21 WRKYs were up-regulated (Table 2). Considering the large number of genes in this family that were up-regulated in our experiment and their importance in the response to pathogen and herbivore attack, we investigated this TF family in more detail. A phylogenetic analysis including WRKY genes from cotton and Arabidopsis was performed (Figure 3 and Additional file 7). The phylogenetic analysis classified the cotton genes in the WRKY clades and identified several putative cotton homologues to Arabidopsis WRKYs that have been previously characterised. Among the genes with an altered response to cotton boll weevil larvae feeding, we identified nine WRKY genes in cotton with putative homologues in Arabidopsis that belong to group III: WRKY46, WRKY30, WRKY64 and WRKY70 (Figure 3 and Additional file 7) [33, 34]. In addition, GhWRKY40-like1 (AT1G80840), GhWRKY33-like1 (AT2G38470) and GhWRKY22-like1 (AT4G01250) genes belong to group IIa, I, and IIe, respectively (Figure 3 and Additional file 7). Other zinc finger-containing transcription factors from the C2H2 and C3H families were also up-regulated during boll weevil larvae feeding. Members of these families are also up-regulated by aphid and Pseudomonas syringae attack in Arabidopsis[35] (Table 2). MYB proteins are central regulators of development, metabolism, and the response to abiotic and biotic stresses [36]. In our analysis, we found six up-regulated genes coding for R2R3-MYB TFs (Table 2). Defence responses regulated by MYB transcription factors promote HR-related cell death and resistance against bacterial and necrotrophic pathogens. MYB transcription factors also play roles in the defence response against insects. Small Cabbage White caterpillars (Pieris rapae) induce local expression of AtMYB102[37]. NAC transcription factors have been shown to regulate plant development, phytohormone signalling and abiotic stress responses [38]. Six NAC genes were differentially expressed in response to cotton boll weevil larvae feeding (Table 2). Some NAC genes are up-regulated in response to feeding by diamond back moths, wounding, and bacterial infection [39, 40]. bHLH transcription factors regulate plant cell and tissue development as well as phytohormone signalling. A limited number of bHLH transcription factors characterised to date have been found to be involved in the defence against pathogens or pests. Nine bHLH TFs were consistently up-regulated in this study. Other differentially expressed (up-regulated) transcription factors belonging to the GRAS family play diverse roles in root and shoot development, gibberellic acid (GA) signalling, phytochrome A signal transduction and the chitin induced defence response [41, 42]. We focused on three members of the AP2/ERE family and nine members of the WRKY family to further characterise the response to cotton boll weevil feeding by qPCR analysis (Table 2).

Figure 3
figure 3

Phylogenetic tree of WRKY domains between cotton and Arabidopsis . The amino acid sequences of the WRKY domain of cotton and Arabidopsis were aligned with MUSCLE, and the phylogenetic tree was constructed using the JTT model with an estimated γ-distribution parameter (G). The maximum-likelihood analyses were performed with the program PhyML version 3.0, and assessment of node confidence was performed using 1,000 bootstrap replicates. The members of group I, II (a-e) and III are labelled according to the classifications of AtWRKY domains by Eulgem et al. [33]. The triangles indicate cotton WRKY genes, and filled triangles represent the genes analysed by qPCR. Contig9787 was named GhWRKY70-like1; was named contig5359 GhWRKY64-like1; and contig16334 was named GhWRKY72-like1 after calculating the p-distance to determine the closest relationship with Arabidopsis members.

Oxidative stress-related genes affected by cotton boll weevil feeding

Oxidative radicals play an important role in plants during various stresses, including biotic stress such as insect infestation [43]. Among the genes up-regulated by oxidative stress, a transcript encoding the respiratory burst oxidative homolog D (RbohD) protein, a key enzyme involved in the reduction of reactive oxygen species (ROS) during the defence response, was identified. In addition, GRX480, a member of the glutaredoxin family that catalyses thiol disulphide reduction and might be a candidate controlling the redox state of regulatory proteins, was also up-regulated. GRX480 might represent a potential regulatory component of SA/JA antagonism [44].

Phytohormone-related genes induced by the boll weevil larvae feeding

Phytohormones are crucial in the stimulation of defence to insect herbivory. Several signalling pathways, including JA, SA, and ET, are believed to orchestrate the induction of insect defences [10]. Genes associated with the biosynthesis, signalling and response to stimuli by JA, SA, ABA (abscisic acid) and ET phytohormones were significantly overrepresented in the flower bud transcriptome affected by the boll weevil based on hypergeometric distribution analysis (p-values ≤ 0.005) (Figure 1). For instance, the genes associated with jasmonate biosynthesis are well represented among the genes up-regulated by cotton boll weevil feeding, including lipoxygenase (LOX), allene oxide synthesis (AOS), and jasmonic acid carboxyl methyltransferase (JMT). We identified genes associated with salicylic acid signalling, a widely known phytohormone involved in the induction of pathogen resistance (PR) genes as well as in the establishment of long-term immunity via the SAR pathways [45] (Figure 1). In addition, other genes up-regulated by boll weevil feeding are involved in ABA catabolism and ET biosynthesis. These include the gene encoding the ABA 8′-hydroxylase (CYP707A3), which catalyses the first step in the oxidative inactivation of ABA [46]. The genes involved in ethylene biosynthesis ACS (1-aminocyclopropane-1-carboxylate synthase) and ACO (1-aminocyclopropane-1-carboxylate oxidase) were also up-regulated after herbivory (Table 3).

Down-regulated genes involved in the response to cotton boll weevil feeding

Several putative genes involved in cell wall formation and degradation were identified as down-regulated. Contig12858 was one of the most highly down-regulated genes in our study (log2-fold change = −34.07, Additional file 2). Contig12858 is similar to genes belonging to the invertase/pectin methylesterase inhibitor family, from which some genes are known to be involved in cell wall modification but have also been reported to play an important role in basal disease resistance [47, 48]. Other genes with potential roles in cell wall reorganisation were down-regulated, such as an endoxyloglucan transferase (contig9845), which has a role in primary cell wall restructuring, as well as an alpha/beta hydrolase (contig12015) and a beta-expansin (contig9652) (Additional file 2). Interestingly, many down-regulated genes encoding heat shock proteins (Hsp) and one heat shock transcription factor (HSF) were observed. Among them, five belong to the small Hsp family, four to the Hsp70 family and six to the Hsp90 family (Additional file 2).

Validation of transcriptome data and evaluation of the spatial expression pattern of defence-related genes in response to infestation by cotton boll weevil larvae using “in situ” qPCR

To better characterise the potential role of some genes in plant defence against insect herbivory, we selected 32 DEGs that responded to feeding by the cotton boll weevil larvae. These genes function in different levels of the biotic stress defence response signalling pathway (Figure 2) and were annotated in biological processes such as the response to chitin and/or cell death (Table 3). Initially, these were analysed with real-time qPCR on 6 mm cotton flower buds infected by larvae or the control samples for further validation of the transcriptome data. The expression patterns observed between the RNA sequencing (RNAseq) and qPCR methods were highly similar (Additional file 8), indicating the accuracy of our transcriptome profile. To verify this observation, we calculated the Pearson correlation coefficient (r) between the different methods for all transcripts, and a correlation coefficient of 0.9406 (P ≤ 0.05) was observed.

To evaluate the spatial expression pattern of these genes by “in situ” qPCR, two different areas of the 6 mm cotton flower buds were isolated by laser microdissection (LMD): the stamens (near the feeding site) and the carpels (away from the injured area) (Figure 4). Among the 28 up-regulated genes selected for analysis, 27 were up-regulated in both organs. Only CML38 showed a contrasting expression pattern of up-regulation in the stamens (near the feeding site) but down-regulation in the carpels (away from the feeding site) (Figure 5). The spatial expression analysis of four down-regulated genes showed that all four were repressed in stamens and carpels (Figure 5). We also evaluated if the expression near the feeding site (stamens) was significantly different from the expression away from the damaged tissue (carpels). Among the 32 genes tested by qPCR, transcripts from 24 genes were more highly expressed near the feeding site (Figure 5). A few genes, such as the transcription factors GhWRKY19-like1, GhWRKY22-like1, GhWRKY40-like1 and ERF-98 and the RbohH and HSP90.1 genes were more highly expressed in carpel tissues compared to stamen tissues (Figure 5 and Additional file 9).

Figure 4
figure 4

Isolation of cotton tissues from paraffin-embedded sections by laser microdissection (LMD). Isolation of cotton tissues from paraffin-embedded sections by laser microdissection (LMD). Sections before LMD (a, c, e), and sections after LM (b, d, f). The area selected for laser microdissection is outlined in green (a region near the damage caused by larvae feeding, which comprised the stamen tissue, viewed at a, c and d) or blue (a region farther from the injured area, which comprised the carpel tissue, viewed at a, e, and f). The assessment of extracted RNA integrity from the stamen and carpel are shown in g and h, respectively. Electropherograms were obtained with an Agilent 2100 Bioanalyser. Open and closed arrowheads indicate the 18S and 28S ribosomal RNA peaks, respectively. RNA quality is expressed as the RNA integrity number (RIN). Scale bars = 100 μm.

Figure 5
figure 5

Comparison of expression levels by “ in situ ” qPCR of a subset of 32 DEG. These genes were examined from two different areas (stamen and carpel) isolated from 6 mm cotton flower buds infested by cotton boll weevil larvae by laser microdissection (LMD) in relation to the control. The reference genes GhACT4 and GhFBX6 were used to normalise the qPCR data. The relative expression level was calculated using the relative expression software tool (REST©), and a subsequent statistical test of the analysed CP values by a Pair-Wise Fixed Reallocation Randomization Test was performed.

The differentially expressed transcriptome of cotton in response to boll weevil larvae feeding

To further evaluate whether cotton and Arabidopsis share a common gene set in response to tissue-chewing pests, we compared the differentially expressed transcriptome with the publicly available microarray data set GEO: GSE10681 (Arabidopsis leaf transcriptome in response to diamond back moth (DBM, Plutella xylostella larvae feeding). We also used a GSEA approach to determine the genes and pathways that were relevant in response to 24 hours of continuous DBM feeding and identified the biological processes that were overrepresented in the ordered list of DEG in this experiment (data not shown). We found several similarities between the cotton transcriptome and the public domain microarray data in terms of enriched biological processes. Biological processes related to hormone biosynthesis and signalling, such as JA, SA, ET and ABA, the responses to organic substances, the negative regulation of programmed cell death, the detection of biotic stimuli, systemic acquired resistance and the regulation of the immune response were overrepresented in both GSEA analyses. However, there were obvious differences in the enriched biological processes, including the response to auxin, gibberellin stimulus and secondary metabolite biosynthetic processes. These biological processes are overrepresented only in the leaf transcriptome analysis of Arabidopsis infested by DMB. On the other hand, the response to chitin and signalling are biological processes overrepresented in cotton flower buds damaged by cotton boll weevil larvae. There were only nine up-regulated genes that overlapped between the cotton larvae in floral tissue and the leaf injured by diamond back moth (Figure 6). Among them, there were transcripts encoding proteins that repress JA signalling, such as JAZ1 (At1g19180), JAZ7 (At3g34600), and JAS1 (At5g13220) as well as a disease resistance protein LRR (At2g34930) that specifically recognises pathogen effectors, resulting in effector-triggered immunity (ETI). Another gene represented in both data sets is a lipoxygenase (LoX4, At1g72520) that dioxygenates unsaturated fatty acids, leading to the lipoperoxidation of biological membranes. LoXs are involved in the apoptosis (programmed cell death) pathway and in biotic and abiotic stress responses in plants [49]. Finally, the gene 2-oxoglutarate and Fe(II)-dependent oxygenase (At5g05600), which is involved in the phenylpropanoid pathway and implicated in the synthesis of defence compounds, was found in both data sets and was also induced by Pieris brassicae eggs and Pieris rapae herbivory [17]. We found only one gene that was down-regulated in both experiments: a germin-like protein (GER3) with oxalate oxidase activity leading to hydrogen peroxide (H2O2) production, which is a ROS that is implicated in the herbivory-induced response in plants [17].

Figure 6
figure 6

Venn diagram comparing the cotton boll weevil ( Anthonomus grandis ) induced transcriptome with the response to another tissue-chewing pest, Plutella xylostella . The number of induced (up-regulated, left column) and repressed (down-regulated, right column) genes after 48 h of cotton boll weevil feeding (A. grandis) was compared to the response of another 24 h herbivory-treatment previously published [39]. Selected genes have a p-value ≤ 0.05 and |logFC| ≥ 2.0.


G. hirsutum L. is an allotetraploid species with a large and complex genome, and comprehensive sequence information describing the genome is still incomplete. Moreover, previous studies have provided limited information on the transcriptional dynamics during the cotton defence response [50]. The recent availability of high-throughput sequencing technologies provides an unprecedented opportunity to thoroughly explore the defence response using large-scale expression profile analysis despite uncharacterised genomic sequences. In the present study, using massive parallel mRNAseq, we report the transcriptional changes in G. hirsutum L. flower buds in response to cotton boll weevil (A. grandis) larvae feeding. To this end, eggs from A. grandis that were ready to hatch were placed on stamens inside 6 mm flower buds, the stage at which the round unicellular microspores are found in the locules. Cotton boll weevil larvae chew for ~48 h on flower tissues, representing the early steps of insect damage. It has been suggested that pollen grains are the initial target of the larvae after egg hatching [51]. Subsequently, the larvae migrate to the ovary where they feed in the ovules and conclude their life cycle, leading to the abscission of floral structures [51].

In this study, we generated over 327 million sequence reads (100 bp in length) of raw sequence data using the Illumina sequencing of two biological replicates from 6 mm cotton flower buds damaged by cotton boll weevil larvae and undamaged flower buds. At the current level of sampling (at least 74 million tags per sample), the sensitivity of the RNASeq approach allowed the accurate identification of differentially expressed genes.

A previous study in cotton showed by RNA-Seq that the transcriptome changes after pathogen (Verticillium dahliae) inoculation, although they generated only 27 million tags of 21 bp in length [50]. Dubey et al. [4] analysed the leaf transcriptome of cotton plants infested by aphids and whiteflies, generating a total of 3.8 million high quality reads through RNA-seq. In our study, the number of up-regulated genes (402) was much greater than the number of down-regulated genes (41) after damage by cotton boll weevil larvae. These results are supported by previous reports showing that larval feeding stimulates induction rather than suppression of genes [18, 39]. Several studies have examined the transcriptional profiles of different modalities of attack, such as pathogen or insect [15, 18, 31]. However, our results contrast with results observed in cotton plants damaged by whiteflies and aphids, in which the number of down-regulated genes was greater than the number of up-regulated genes in both cases after damage [4].The MapMan analysis of the DEGs generated a representative overview of the plant response after larvae damage (Figure 2). This analysis suggested that the cotton plant could perceive boll weevil larvae-derived physical and chemical cues such as compounds in the oviposition fluid, chitin, and insect oral secretions. These elicitors dramatically alter the expression of several genes in the plant immune response pathways that are involved in biotic stress response. The herbivory-induced changes are mediated by elaborate signalling networks, including receptors/sensors, Ca2+ influx, kinase cascades, transcription factors, reactive oxygen species and phytohormone signalling pathways. We were able to identify genes belonging to all major networks previously described and thus relate these to the plant response to herbivory.

The plant immune system is activated by various herbivore (or microbial)-associated molecular patterns (HAMP or MAMP) that are detected as non-self molecules. Such patterns are recognised by immune receptors that are either cytoplasmic or localised on the plasma membrane. Cell surface receptors include receptor-like kinases (RLK) that frequently contain extracellular leucine-rich repeats and an intracellular kinase domain for the activation of downstream signalling, as well as receptor-like proteins (RLP) that lack this signalling domain. It is therefore hypothesised that RLKs are required for RLPs to activate downstream signalling [52]. Several putative receptor-like kinases, such as the homologues of CCR4, RIPK and EVR, were up-regulated after cotton boll weevil larvae damage (Table 3). CCR4 encodes a transmembrane protein kinase (RLK) previously characterised as a key gene in the plant-insect interaction. Little et al. [17] showed that 41 RLKs, including CCR4, were induced by oviposition by P. brassicae, suggesting that RLKs play an important role in the detection of herbivore-associated molecular patterns. Interestingly, CCR4 was significantly more highly expressed in the stamen close to the larval damage than in the carpel, as assessed by our in situ qPCR analysis (Figure 5). Furthermore, the RIPK gene encodes a cytoplasmic NB-LRR immune receptor that recognises the AvrB and AvrRpm1 effectors from Pseudomonas syringae, leading to the activation of ETI [53]. This gene was also more highly expressed in the stamen than in the carpel tissue in our in situ qPCR analysis. Another leucine rich repeat transmembrane protein kinases analysed in our study, EVR, was equally induced by boll weevil larvae in the stamen and the carpel, is also expressed in response to P. syringae, and is involved in the regulation of cell death and innate immunity (Figure 5) [52].

The recognition of conserved HAMP, such as chitin or fatty acid-amino acid conjugates (FAC, the major component of insect oral secretions) by PRR triggers intracellular signalling via a mitogen-activated protein kinase (MAPK, MAPKK, MAPKKK) cascade [27]. Our study identified homologues to MAPK3 (At3g45640) and MKK9 (At1g73500) genes with the GO term “response to chitin”. The in situ qPCR analysis showed that these genes were highly expressed in the stamen and carpel from cotton flower buds infested by boll weevil larvae (Figure 5). MAPK, MAPKK, and MAPKKK activate several transcription factors such as ERFs and WRKYs during the response to pest attacks, as has been shown for ERF-5[54]. A previous study in Arabidopsis showed that ERF-5 is phosphorylated by MAPK3 (At3g45640), suggesting that ERF-5 might play an important role in plant defence by positively regulating ethylene (ET) signalling [54]. The GhERF-4 and GhERF-5 genes were more highly expressed in the stamen than in the carpel, and the GhERF-98 gene was slightly more highly expressed in the carpel (Figure 5). ERF-4, ERF-5 and ERF-98 belong to the AP2/EREB family and were previously identified in the response to nematode or fungal inoculation (Alternaria brassicicola) as well as in response to chitin, a main elicitor of the plant defence response against insects [31, 55].

RNAseq analysis identified transcripts for 88 TFs that were differentially expressed in 6 mm cotton flower buds damaged by boll weevil larvae. In the present study, we selected the 21 cotton genes identified as WRKY for a phylogenetic analysis in combination with the full set of Arabidopsis WRKY genes. The expression patterns of nine WRKY genes were explored in detail with a combination of LMD and qPCR. The transcript level of GhWRKY30-like1 gene was the same in stamen and carpel tissue and was induced relative to the control. GhWRKY33-like, GhWRKY46-like1, GhWRKY64-like1, GhWRKY70-like1 and GhWRKY72-like1 were more highly expressed near the damage in stamen tissues. In Arabidopsis, it was also shown that WRKY46 is specifically induced by salicylic acid (SA) and infection by the biotrophic pathogen P. syringae, working redundantly with the structurally related genes WRKY70 and WRKY53 to positively regulate basal resistance to P. syringae[56]. Previous results have shown that the Arabidopsis WRKY genes are rapidly and strongly induced by chitin [31]. Some WRKYs are target genes of the upstream MAPK3 (At3g45640), such as WRKY33. Studies in Arabidopsis have shown that in response to pathogen invasion, MAPK3 (At3g45640) phosphorylates WRKY33, which directly binds to the W-boxes in the promoter of the ACS6 gene in vivo, suggesting that WRKY33 is directly involved in the activation of ACS6 expression downstream of the MAPK3 cascade [57]. This signalling event produces high levels of ethylene, which plays an important role in plant immunity. Moreover, AtWRKY33 is required to activate the synthesis of antimicrobial substances such as camalexin [34]. Interestingly, the other three WRKY genes (GhWRKY19-like1, GhWRKY22-like1 and GhWRKY40-like1) showed a higher expression in the carpel than in the stamen. In an elegant series of experiments, Asai et al. [58] showed that AtWRKY22 functions downstream of the flagellin receptor, a leucine-rich repeat (LRR) receptor kinase, and that a mitogen-activated protein kinase (MAPK) cascade conferred resistance to bacterial (P. syringae) and fungal (Botrytis cinerea) pathogens. These results strongly indicate that signalling events initiated by boll weevil larvae feeding might converge on a conserved system of pathogen-associated molecular pattern recognition, a MAPK cascade and transcription factors.

Relatively little is known about the signal transduction pathways triggered by insect damage. Calcium ions (Ca2+) have been implicated as a second messenger in many plant signalling pathways including responses to herbivory [12]. It has been shown that mechanical wounding and insect damage trigger a transient increase in cytosolic Ca2+ levels [59]. For instance, the Egyptian cotton worm (Spodoptera littoralis) caused a transient increase in cytosolic Ca2+ in Phaseolus lunatus cells adjacent to the insect damage [59]. Cotton genes encoding calmodulin-like proteins (EDA39, CML38 and CML24) showed a higher expression close to the damage (stamen tissues) than in carpels, suggesting that these proteins might be a component of Ca2+ signalling that modulate plant defence responses against herbivory attack.

Several studies have revealed that ROS are implicated in herbivory-induced responses in plants [12, 43, 60]. Superoxide anion (O2−), hydrogen peroxide (H2O2), singlet oxygen (1O2), and hydroxyl radical (·OH) are collectively called ROS; they are produced in mitochondria, chloroplasts, and peroxisomes, as well as on the external surfaces of plasma membranes. Feeding of Helicoverpa zea on soybean, as well as S. littoralis and Tetranychus urticae feeding on Medicago truncatula plants, increase ROS levels [12, 61]. Our analysis identified several transcripts encoding enzymes directly involved in oxidative stress that are up-regulated in response to feeding by boll weevil larvae. In situ qPCR analysis showed that transcripts encoding an NADPH oxidase (RbohD) were expressed at a higher level in carpel tissues than in stamen tissues compared to the control samples (Figure 5). AtrbohD was largely responsible for the accumulation of ROS during the defence response in Arabidopsis against P. syringae and Peronospora parasitica[60]. Our results suggest that the oxidative burst is also part of the cotton plant’s response to boll weevil larvae feeding.

The boll weevil larval damage in flower buds also induced the expression of genes associated with phytohormone biosynthetic processes and signalling. Among the gene transcripts associated with JA, ET and ABA biosynthesis selected for analysis by qPCR were AOS, ACS, ACO and CYP701A3. All of these genes showed a significantly higher expression in the stamen and carpel tissue compared to the control samples. Among the JA metabolism genes that modulate plant responses to biotic stress by boll weevil larvae are AOS, which encodes the second enzyme involved in JA biosynthesis, and JMT, which is involved in the synthesis of the volatile compound methyl-JA (MeJA) used for plant defence. Jasmonate plays a central role in regulating the defence responses to herbivores that inflict various types of tissue damage. Jasmonate mutants are affected in their resistance to a wide range of arthropod herbivores, including caterpillars (Lepidoptera), beetles (Coleoptera), thrips (Thysanoptera) and leafhoppers (Homoptera) [9, 59, 62]. Jasmonate has been associated with the synthesis of proteinase inhibitors (PI), defence-related volatile compounds and secondary metabolites, such as nicotine, active phenolics and phytoalexins [12]. The genes involved in the ET biosynthesis pathway, ACS and ACO, were also up-regulated by boll weevil larvae damage. The primary function of ET in plant resistance to herbivores is the fine-tuning of JA-induced responses [12]. In tomato plants, ET potentiates the JA-induced transcript accumulation of secondary metabolites such as PIs [63]. The treatment of Arabidopsis plants with ethephon, a synthetic ET donor, transiently elevates the levels of JA and AOS transcripts [64]. Blocking the perception of ET with 1-MCP diminishes herbivory-induced volatile emission, which is mainly regulated by JA [65]. By using plants that ectopically express a loss-of-function ETR gene, von Dahl et al. [66] showed that ET signalling is crucial for increasing basal levels of nicotine, an effective defence against herbivory. These results indicated that JA and ET may be involved in elevating the basal resistance of cotton plants to herbivory.

The ubiquitin-ligase genes are also known to play a role in the plant defence response [67]. Yang et al. [68] showed that the E3 ligase activity of Arabidopsis PUB17 is required for the initiation of the hypersensitive response (HR). Therefore, two predicted members of this family, PUB23 and CMPG1, which are among the DEGs in our experiment and were annotated in response to chitin and/or death, were also analysed by in situ qPCR. Spatial analysis of gene expression revealed that the transcripts of both genes were strongly up-regulated at the site close to the damage (stamen) compared to carpel tissues. These genes are among the 30 ubiquitin-ligase genes that up-regulated either 15 or 30 min after chitooctaose treatment in Arabidopsis[31]. Previous studies have shown that the expression of the PUB23 gene in Arabidopsis was induced by PAMPs and infection by pathogens. Moreover, the pub22/pub23/pub24 triple mutant displayed negative regulation of PAMP-triggered immunity [69].

Among the repressed gene transcripts related to the herbivory response, we identified transcripts for one germin-like protein (GER3) and many genes encoding heat shock proteins (Hsp). qPCR analysis revealed that the GER3 transcript was strongly down-regulated in the stamen and in the carpel (Figure 5). Interestingly, GER3 is included in the overlap between the down-regulated genes of the leaf transcriptome analysis of Arabidopsis damaged by DMB larvae and our study (Figure 6). GER3 leads to hydrogen peroxide (H2O2) production, a type of ROS that is implicated in the herbivory-induced response in plants, and was also repressed by P. brassicae eggs and P. rapae herbivory in a study in A. thaliana[17]. The transcriptional repression of genes encoding the heat shock proteins HSP20, HSP90.1 and HSP90.2 was observed in both stamen and carpel tissues (Figure 5). In many plant species, HSP90 isoforms are required for disease resistance against invading pathogens. HSP90 interacts with the disease resistance signalling components SGT1b and RAR1 and is required for RPS2-mediated resistance, a disease resistance (R) intracellular protein that specifically recognises pathogen effectors. For example, the AtHSP90.1 and AtHSP90.2 genes in Arabidopsis are required for RPS2-mediated resistance against P. syringae expressing AvrRpt2 and for RPM1-mediated resistance to P. syringae expressing AvrRPM1, respectively [7072]. HSP90 is also essential for Rx-mediated resistance to Potato virus X (PVX), N-mediated resistance to Tobacco mosaic virus, and Pto-mediated resistance to P. syringae expressing AvrPto[73, 74]. In contrast, the hsp90.2-3 mutant with a point mutation in the ATP-binding domain of AtHSP90.2 is known to be more sensitive to biotrophic pathogens. In our results, several hsp genes were down-regulated, suggesting that larvae attempt to overcome the plant immune defence triggered by the effector.


Alternative strategies for protecting crops from insect pests are exploring the endogenous resistance mechanisms exhibited by plants to most herbivore insects through a greater understanding of induced defences in plants. This study aimed to unravel the changes in the transcriptome of G. hirsutum flower buds in response to feeding by cotton boll weevil (A. grandis) larvae. A large number of cotton transcripts and biological processes were significantly altered upon larvae infestation. Among the changes observed, we highlighted the induction of transcripts for the receptors/sensors that recognise chitin, insect oral secretions and signals from injured plant cells; differential modulation of transcripts encoding enzymes related to kinase cascades; transcription factors (such as WRKY and ERF); Ca2+ influx; reactive oxygen species; and modulation of transcripts encoding enzymes from phytohormone signalling pathways, mainly the JA and ET pathways. Cotton boll weevil larvae feeding affected many genes that have also been shown to be regulated in response to microbial or fungal infection, indicating the existence of complex crosstalk in the response to these different pathogens. The qPCR analysis associated with microdissection showed that almost all of the genes were more highly expressed in the stamen, the region close to the damage caused by larvae feeding, than in the carpel, farther away from the injured area. However, herbivory-induced defence responses in cotton flower buds were observed not only in the wounded regions but also in undamaged regions of the damaged flower buds. It is possible that a signal travels to other parts of the flower bud that transmit an herbivory alert. Although the identity of the mobile signal responsible for this systemic response remains unknown in many plants, studies performed in solanaceous plants demonstrated that herbivory rapidly induces a short-distance mobile signal that travels and activates MAPKs and triggers JA accumulation. Deciphering the signals that regulate herbivore-responsive gene expression in G. hirsutum flower buds might provide new strategies to manipulate this response for the production of insect-resistant transgenic plants.


Plant material and A. grandisinfestation assay

Cotton (Gossypium hirsutum L. variety “BRS Cedro”) plants were grown under a controlled temperature (27 ± 2°C) and natural photoperiod at Embrapa Genetic Resources and Biotechnology in Brasilia (DF, Brazil). Three-month old plants containing 6 mm flower buds were selected for the experiments. A population of A. grandis (Coleoptera: Curculionidae) was maintained at 27 ± 2°C, 70 ± 10% relative humidity, and a photoperiod of 14 h. Insects were routinely maintained on a standard rearing diet [75]. A 6 mm cotton flower bud previously drilled with a needle was inoculated with an A. grandis egg containing an active embryo. The orifice resulting from drilling was sealed with Vaseline to prevent dehydration of the egg. The hatching of larvae occurred approximately 2 h after egg inoculation. The larvae were removed with tweezers from the flower buds 48 h after inoculation, and four flower buds from different plants were immediately frozen in liquid nitrogen for the isolation of total RNA to perform RNAseq experiments. Similarly, a set of twelve flower buds infected with larvae were collected and fixed in cold ethanol:acetic acid (3:1) to subsequently isolate the selected tissue through laser microdissection. The insect infection experiments were performed in two biological replicates. The control samples were drilled cotton flower buds into which no eggs were introduced.

RNA isolation and sample preparation

Total RNA extraction was performed from 100 mg of macerated plant tissue in liquid nitrogen using the Invisorb Spin Plant RNA Mini kit (Invitek) according to the manufacturer’s protocol. RNA quality and quantity were determined using a Nanodrop 2000 (Thermo Scientific) and a Bioanalyser Chip RNA 7500 series II (Agilent), respectively. The analysis of RNA integrity was judged by the RNA Integrity Number (RIN), which was calculated with 2100 Expert Software. The software algorithm, which was developed by Schroeder et al. [76], categorises total RNA quality on a scale from 1 to 10, in which 10 corresponds to intact RNA and 1 corresponds to highly degraded RNA. Generally, plant RNA with a RIN value of 6–7 was of acceptable quality for qPCR and RNAseq analyses. RNA Integrity Number (RIN) values were greater than 9.0 for all samples. Two biological replicates of inoculated flower buds and controls were used for transcriptome sequencing. Library preparation and massive parallel sequencing were performed by Eurofins MWG Operon (Huntsville, AL). Sequencing libraries were prepared using NEBNext Ultra Directional RNA Library Prep Kit (New England Biolabs, MA), which included Poly-A-containing mRNA isolation from 5 μg total RNA using two rounds of purification with poly-T oligo-attached magnetic beads and fragmentation by sonication. First strand cDNA was generated using reverse transcriptase and random primers. Following the second strand cDNA synthesis and adaptor ligation, 400 bp cDNA fragments were isolated using gel electrophoresis and amplified by PCR. The products were loaded onto an Illumina HiSeq™ 2000 instrument and subjected to 200 cycles of paired-end (2 × 100 bp) sequencing. The processing of fluorescent images into sequences, base-calling, and quality value calculations were performed using the Illumina data processing pipeline (version 1.8).

Mapping of short reads and assessment of differential gene expression

Raw reads were filtered to obtain high-quality reads by removing low-quality reads containing more than 30% bases with Q < 20. After trimming low-quality bases from the 5′ and 3′ ends of the remaining reads, the resulting high-quality reads were aligned against a G. hirsutum expressed sequence tag (EST) assembly (with 28,432 unique genes/contigs) as the transcriptome reference (cotton EST database; [77] using the BWA package [78]. Differential expression was estimated and tested with the DEGseq software package in R-bioconductor 2.15 for each library with reference to the control [79]. The count data were normalised to the total number of counts, accounting for the variance and the mean of the biological replicates. Contigs with adjusted p-value ≤ 0.05 and estimated absolute log2 (FC) ≥ 2 were determined to be significantly differentially expressed genes (DEG) and were selected.

Functional annotation

The DEG were subjected to blast using the program Blastx with the TAIR 9 protein database and further classified into categories according to the GO classification system (Additional file 2). The Blast2GO program ( [80] was also used to confirm the annotation of transcripts (Additional file 2). RNA-seq data can be accessed at the NCBI bioprojects under the accession number PRJNA245406.

To identify relevant molecular mechanisms potentially associated with the response to cotton boll weevil larvae development within cotton flower buds, GSEA (Gene Set Enrichment Analysis) was performed [81]. A gene set was defined as all DEG, with annotations according to A. thaliana, that share the same ontology based on the GO database. The GSEA method identified biological processes (BP), molecular functions (MF) and cellular components (CC) that were overrepresented among the list of DEG (Figure 1 and Additional file 3). The overrepresentation was assessed with a statistical score based on a hypergeometric test with p-values ≤ 0.005.

In addition, the differentially expressed genes were also functionally analysed using MapMan software, which is a user-driven tool that displays large genomic datasets onto diagrams of metabolic pathways or other processes such as biotic stress [82] (Figure 2, adapted by [10]). Finally, the differentially expressed genes were compared with the public databases generated from microarray analysis of Arabidopsis thaliana leaves that were infested by diamond back moth (Plutella xylostella) larvae 24 h after the onset of herbivory (GEO:GSE10681) (Figure 6) [39].

Defence-related genes in response to infestation by cotton boll weevil larvae used for qPCR analysis

Real time PCR (qPCR) assays were performed on 32 differentially expressed genes identified by RNAseq analysis (Additional file 9). Twenty-eight of them were reported as up-regulated in cotton boll weevil-infested flower buds: three receptor-like kinase proteins (CCR4, RIPK and EVR); two proteins from the mitogen-activated kinase signalling cascade (MAPK3, MKK9); three calmodulin-like proteins (CML28, CML24, EDA39); two genes from the oxidative burst process (GRX480, RbohD); twelve TFs, including nine WRKY transcription factors (WRKY19, 22, 30, 33, 40, 46, 64, 70 and 72) and three AP2/ERE transcription factors (ERF-4, ERF-5, ERF-98); four genes involved in phytohormone biosynthesis (CPY707A3, AOS, ACS, ACO); and two ubiquitin ligase proteins involved in targeted protein degradation (PUB23 and CMPG1). Among the down-regulated DEGs, we evaluated the expression pattern of a germin-like protein (GER3) and three members of the heat shock protein family (HSP90.1, HSP90.2 and HSP20). Primers were designed using Primer 3 software tools ( with a melting temperature of 60°C, an amplicon length of 150 to 200 bp, and a GC content of 50 to 60% [83].

Tissue fixation, embedding, sectioning, and laser microdissection

To assess the spatial expression pattern of the genes described above, tissues from two different areas were isolated from 6 mm cotton flower buds infested by cotton boll weevil larvae using laser microdissection (LMD). These areas (Figure 4) included the following: (a) a region near the damage caused by larvae feeding, consisting of stamen tissues and (b) another region farther away from the injured area, consisting of carpel tissue. The 6 mm cotton flower buds infested by larvae and control samples were fixed overnight at 4°C in Farmer’s Fixative (3:1 ethanol:acetic acid) and dehydrated. After fixation, samples were prepared as previously described [84]. Briefly, the samples were subjected to an ethanol/xylene series followed by a xylene/paraplast series before being embedded in paraplast X-tra (Sigma-Aldrich). Embedded samples were sectioned using a rotatory microtome (HM 325, Microm). The 16-μm-thick sections from the infected flower buds and control samples were mounted on Leica FrameSlides PET-Membrane 1.4 μm slides with methanol (No. 11505151). The slides were dried at 40°C for 5 min and then stored at 4°C until laser microdissection was performed. Laser microdissection was performed using a Leica LDM CTR 6500.

RNA preparation, RNA amplification, and “in situ” qPCR

Isolated areas of carpel and stamen with biological replicates of infested and control samples (total of 180–200 microdissected cuts for each sample) were collected in microfuge tubes containing 100 μl RNALater solution (Ambion) and were used for RNA isolation (Figure 4). RNA was extracted from the captured tissues using the RNeasy Plant Mini kit (Qiagen) according to the manufacturer’s protocol, including a 15 min DNase incubation step (RNase-free DNase Set, Qiagen). The eluted RNA was concentrated under a vacuum until the remaining volume was approximately 10 μl. The quality of total RNA extracted from LMD-collected tissues was assessed using an RNA 6000 Pico kit on the Agilent 2100 Bioanalyser (Agilent Technologies). RNA Integrity Number (RIN) values were greater than 6.0 for all samples. The isolated total RNA was amplified using a MessageAmp aRNA Kit (Ambion) following the manufacturer’s instructions. One round of amplifications was performed, and the resulting anti-sense RNA (aRNA) was quantified by OD260. Equal amounts of aRNA (5.0 μg) of the two samples in biological replicates were reverse transcribed using 0.5 μl of Random Primers (C1181, Promega) and Superscript III reverse transcriptase (Invitrogen) according to the manufacturer’s instructions. Polymerase chain reactions were performed in the 7500 Fast Real-Time PCR detection system (Applied Biosystem) and SYBR® Green was used to monitor dsDNA synthesis. PCR efficiencies and the optimal quantification cycle threshold or Cq values were estimated using the online Real time PCR Miner[85]. Two independent biological samples for each experimental condition were evaluated using triplet technical replicates. The reference genes GhACT4 and GhFBX6 were used to normalise the qPCR data and have been discussed previously [86] (Additional file 9). The relative expression level was calculated using the relative expression software tool (REST©), which compares two treatment groups with multiple data points in the sample compared to the control groups and calculates the relative expression ratio between them. The mathematical model used is based on the PCR efficiencies and the mean CP deviation between the sample and control group of target genes and is normalised to the mean CP deviation of the reference genes [83]. The advantage of REST is the use of a subsequent statistical test of the analysed CP values by a Pair-Wise Fixed Reallocation Randomization Test[87].

Phylogenetic analysis of WRKY genes involved in the response to cotton boll weevil feeding

To identify the twenty-one cotton WRKY genes that were differentially expressed in our RNAseq experiment, a phylogenetic tree based on the WRKY domain of Arabidopsis was conducted (Additional files 2 and 7). The WRKY domain boundary was defined as described in Eulgem et al. [33] (excluding the N-terminal domains of Group I). Motif detection was performed with MEME 4.0 software [88] and the multiple amino acid sequence alignment of WRKY domains were conducted using MUSCLE ( from 72 AtWRKY protein sequences from Arabidopsis (downloaded from AtTFDB ( and 21 cotton WRKY genes. Phylogenetic trees were constructed after using a model test to pick the most appropriate evolutionary model for our analysis in the Mega 5.2 program [89]. The best-fitting amino acid substitution model for these gene families was the JTT model [90] with an estimated γ-distribution parameter (G). The maximum-likelihood analyses were performed with the program PhyML version 3.0 ( [91] and assessment of node confidence was performed using 1,000 bootstrap replicates. The trees were visualised and optimised in the Mega 5.2 program (Figure 3).


  1. Wendel JF, Cronn RC: Polyploidy and the evolutionary history of cotton. Adv Agron. 2003, 78: 139-186.

    Article  Google Scholar 

  2. Lee JJ, Woodward AW, Chen ZJ: Gene expression changes and early events in cotton fibre development. Ann Bot. 2007, 100 (7): 1391-1401. 10.1093/aob/mcm232.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  3. Razaq M, Aslam M, Shad SA, Naeem M: Evaluation of some new promising cotton strains against bollworm complex. Evaluation. 2004, 15 (3): 313-318.

    Google Scholar 

  4. Dubey NK, Goel R, Ranjan A, Idris A, Singh SK, Bag SK, Chandrashekar K, Pandey KD, Singh PK, Sawant SV: Comparative transcriptome analysis of Gossypium hirsutum L. in response to sap sucking insects: aphid and whitefly. BMC Genomics. 2013, 14: 241-261. 10.1186/1471-2164-14-241.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  5. Greenberg SM, Sappington TW, Setamou M, Coleman RJ: Influence of different cotton fruit sizes on boll weevil (Coleoptera: Curculionidae) oviposition and survival to adulthood. Environ Entomol. 2003, 33: 443-449.

    Article  Google Scholar 

  6. Martins WFS, Ayres CFJ, Lucena WA: Genetic diversity of Brazilian naturalpopulations of Anthonomus grandis Boheman(Coleoptera: Curculionidae), the major cottonpest in the New World. Genet Mol Res. 2007, 6: 23-32.

    CAS  PubMed  Google Scholar 

  7. Oliveira-Neto OB, Batista JAN, Rigden DJ, Franco OL, Fragoso RR, Monteiro ACS, Monnerat RG, Grossi-de-Sa MF: Molecular cloning of a cysteine proteinase cDNA from the cotton boll weevil Anthonomus grandis (Coleoptera: Curculionidae). 2004 Jun;68(6):1235–42. Biosci Biotechnol Biochem. 2004, 68 (6): 1235-1242. 10.1271/bbb.68.1235.

    Article  PubMed  Google Scholar 

  8. Dicke M, van Poecke RMP, Boer JG: Inducible indirect defence of plants: from mechanisms to ecological functions. Basic Appl Ecol. 2003, 4: 27-42. 10.1078/1439-1791-00131.

    Article  CAS  Google Scholar 

  9. Kessler A, Baldwin IT: Plant responses to insect herbivory: the emerging molecular analysis. Annu Rev Plant Biol. 2002, 53: 299-328. 10.1146/annurev.arplant.53.100301.135207.

    Article  CAS  PubMed  Google Scholar 

  10. Erb M, Meldau S, Howe GA: Role of phytohormones in insect-specific plant reactions. Trends Plant Sci. 2012, 17: 250-259. 10.1016/j.tplants.2012.01.003.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  11. Jones JDG, Dangl JL: The plant immune system. Nature. 2006, 444: 323-329. 10.1038/nature05286.

    Article  CAS  PubMed  Google Scholar 

  12. Wu J, Baldwin IT: New insights into plant response to the attack from insect herbivore. Annu Rev Genet. 2010, 44: 1-24. 10.1146/annurev-genet-102209-163500.

    Article  CAS  PubMed  Google Scholar 

  13. Ferry N, Edwards MG, Gatehouse JA, Gatehouse AMR: Plant–insect interactions: molecular approaches to insect resistance. Curr Opin Biotechnol. 2004, 15: 155-161. 10.1016/j.copbio.2004.01.008.

    Article  CAS  PubMed  Google Scholar 

  14. Wei X, Zhang X, Shen D, Wang H, Wu Q, Lu P, Qiu Y, Song J, Zhang Y, Li X: Transcriptome analysis of Barbarea vulgaris infested with diamondback moth (Plutella xylostella) larvae. PlosOne. 2013, 8: 1-19.

    Google Scholar 

  15. Mafra V, Martins PK, Francisco CS, Ribeiro-Alves M, Freitas-Astúa J, Machado MA: Candidatus Liberibacter americanus induces significant reprogramming of the transcriptome of the susceptible citrus genotype. BMC Genomics. 2013, 14: 247-10.1186/1471-2164-14-247.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  16. Xu L, Zhu L, Tu L, Liu L, Yuan D, Jin L, Long L, Zhang X: Lignin metabolism has a central role in the resistance of cotton to the wilt fungus Verticillium dahliae as revealed by RNA-Seq-dependent transcriptional analysis and histochemistry. J Exp Bot. 2011, 62 (15): 5607-5621. 10.1093/jxb/err245.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  17. Little D, Gouhier-Darimont C, Bruessow F, Reymond P: Oviposition by pierid butterflies triggers defense responses in Arabidopsis. Plant Physiol. 2007, 143: 784-800.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  18. De Vos M, Oosten VRV, van Poecke RMP, Van Pelt JA, Pozo MJ, Mueller MJ, Buchala AJ, Métraux JP, Van Loon LC, Dicke M, Pieterse CMJ: Signal signature and transcriptome changes of Arabidopsis during pathogen and insect attack. MPMI. 2005, 18: 923-937. 10.1094/MPMI-18-0923.

    Article  CAS  PubMed  Google Scholar 

  19. Garg R, Patel RK, Tyagi AK, Jain M: De novo assembly of chickpea transcriptome using short reads for gene discovery and marker identification. DNA Res. 2011, 18: 53-63. 10.1093/dnares/dsq028.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  20. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Baren MJ, Salzberg SL, Wold BJ, Pachter L: Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010, 28: 511-515. 10.1038/nbt.1621.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  21. Autran D, Baroux C, Raissig MT, Lenormand T, Wittig M, Grob S, Steimer A, Barann M, Klostermeier UC, Leblanc O, Vielle-Calzada JP, Rosenstiel P, Grimanelli D, Grossniklaus U: Maternal epigenetic pathways control parental contributions to Arabidopsis early embryogenesis. Cell. 2011, 145: 707-719. 10.1016/j.cell.2011.04.014.

    Article  CAS  PubMed  Google Scholar 

  22. Morin RD, O’Connor MD, Griffith M, Kuchenbauer F, Delaney A, Prabhu AL, Zhao Y, McDonald H, Zeng T, Hirst M, Eaves CJ, Marra MA: Application of massively parallel sequencing to microRNA profiling and discovery in human embryonic stem cells. Genome Res. 2008, 18: 610-621. 10.1101/gr.7179508.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  23. Tarazona S, García-Alcalde F, Dopazo J, Ferrer A, Conesa A: Differential expression in RNA-seq: A matter of depth. Genome Res. 2011, 21: 2213-2223. 10.1101/gr.124321.111.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  24. Garber M, Grabherr MG, Guttman M, Trapnell C: Computational methods for transcriptome annotation and quantification using RNA-seq. Nat Methods. 2011, 8: 469-477. 10.1038/nmeth.1613.

    Article  CAS  PubMed  Google Scholar 

  25. Pavithra SR, Kumar R, Tatu U: Systems analysis of chaperone networks in the malarial parasite Plasmodium falciparum. PLoS Comput Biol. 2007, 14: 1701-1715.

    Google Scholar 

  26. Panstruga R, Parker EJ, Schulze-Lefert P: SnapShot: plant immune response pathways. Cell. 2009, 6: 978-

    Google Scholar 

  27. He P, Shan L, Lin NC, Martin GB, Kemmerling B, Nurnberger T, Sheen J: Specific bacterial suppressors of MAMP signaling upstream of MAPKKK in arabidopsis innate immunity. Cell. 2006, 125: 563-575. 10.1016/j.cell.2006.02.047.

    Article  CAS  PubMed  Google Scholar 

  28. Kannan P, Pandey D, Gupta AK, Punetha H, Taj G, Kumar A: Expression analysis of MAP2K9 and MAPK6 during pathogenesis of Alternaria blight in Arabidopsis thaliana ecotype Columbia. Mol Biol Rep. 2012, 39: 4439-4444. 10.1007/s11033-011-1232-1.

    Article  CAS  PubMed  Google Scholar 

  29. Takahashi F, Mizoguchi T, Yoshida R, Ichimura K, Shinozaki K: Calmodulin-dependent activation of MAP kinase for ROS homeostasis in arabidopsis. Mol Cell. 2011, 41: 649-660. 10.1016/j.molcel.2011.02.029.

    Article  CAS  PubMed  Google Scholar 

  30. Ma W, Smigel A, Tsai Y, Braam J, Berkowitz GA: Innate immunity signaling: cytosolic Ca2+ elevation is linked to downstream nitric oxide generation through the action of calmodulin or a calmodulin-like protein. Plant Physiol. 2008, 148: 818-828. 10.1104/pp.108.125104.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  31. Libault M, Wan J, Czechowski T, Udvardi M, Stacey G: Identification of 118 Arabidopsis transcription factor and 30 ubiquitin-ligase genes responding to chitin, a plant-defense elicitor. Mol Plant Microbe Interact. 2007, 20: 900-911. 10.1094/MPMI-20-8-0900.

    Article  CAS  PubMed  Google Scholar 

  32. Rushton PJ, Somssich IE, Ringler P, Shen QJ: WRKY transcription factors. Trends Plant Sci. 2010, 15: 247-258. 10.1016/j.tplants.2010.02.006.

    Article  CAS  PubMed  Google Scholar 

  33. Eulgem T, Somssich IE: Networks of WRKY transcription factors in defense signaling. Curr Opin Plant Biol. 2007, 10: 366-371. 10.1016/j.pbi.2007.04.020.

    Article  CAS  PubMed  Google Scholar 

  34. Pandey SP, Somssich IE: The role of WRKY transcription factors in plant immunity. Plant Physiol. 2009, 150: 1648-1655. 10.1104/pp.109.138990.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  35. Barah P, Winge P, Kusnierczyk A, Tran DH, Bones AM: Molecular signatures in arabidopsis thaliana in response to insect attack and bacterial infection. Plos One. 2013, 8: 1-24.

    Article  Google Scholar 

  36. Kaur H, Heinzel N, Schottner M, Baldwin IT, Galis I: R2R3-NaMYB8 regulates the accumulation of phenylpropanoid-polyamine conjugates, which are essential for local and systemic defense against insect herbivores in Nicotiana attenuata. Plant Physiol. 2010, 152: 1731-1747. 10.1104/pp.109.151738.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  37. Verk MC, Gatz C, Linthorst HJM: Transcriptional regulation of plant defense responses. Adv Bot Res. 2009, 51: 397-438.

    Article  Google Scholar 

  38. Nakashima K, Tran LSP, Nguyen DV, Fujita M, Maruyama K, Todaka D, Ito Y, Hayashi N, Shinozaki K, Yamaguchi-Shinozaki K: Functional analysis of a NAC-type transcription factor OsNAC6 involved in abiotic and biotic stress-responsive gene expression in rice. Plant J. 2007, 51: 617-630. 10.1111/j.1365-313X.2007.03168.x.

    Article  CAS  PubMed  Google Scholar 

  39. Ehlting J, Chowrira SG, Mattheus N, Aeschliman DA, Arimura GI, Bohlmann J: Comparative transcriptome analysis of Arabidopsis thaliana infested by diamondback moth (Plutella xylostella) larvae reveals signatures of stress response, secondary metabolism, and signalling. BMC Genomics. 2008, 9: 154-10.1186/1471-2164-9-154.

    Article  PubMed Central  PubMed  Google Scholar 

  40. Collinge M, Boller T: Differential induction of two potato genes, Stprx2 and StNAC, in response to infection by Phytophthora infestans and to wounding. Plant Mol Biol. 2001, 46: 521-529. 10.1023/A:1010639225091.

    Article  CAS  PubMed  Google Scholar 

  41. Hirsch S, Oldroyd GED: GRAS-domain transcription factors that regulate plant development. Plant Signal Behav. 2009, 8: 698-700.

    Article  Google Scholar 

  42. Day RB, Shibuya N, Minami E: Identification and characterization of two new members of the GRAS gene family in rice responsive to N-acetylchitooligosaccharide elicitor. Biochim Biophys Acta. 2003, 1625: 261-268. 10.1016/S0167-4781(02)00626-7.

    Article  CAS  PubMed  Google Scholar 

  43. Liu X, Williams CE, Nemacheck JA, Wang H, Subramanyam S, Zheng C, Chen MS: Reactive oxygen species are involved in plant defense against a gall midge. Plant Physiol. 2010, 152: 985-999. 10.1104/pp.109.150656.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  44. Ndamukong I, Abdallat AA, Thurow C, Fode B, Zander M, Weigel R, Gatz C: SA-inducible Arabidopsis glutaredoxin interacts with TGA factors and suppresses JA-responsive PDF1.2 transcription. Plant J. 2007, 50: 128-139. 10.1111/j.1365-313X.2007.03039.x.

    Article  CAS  PubMed  Google Scholar 

  45. Smith JL, De Moraes CM, Mescher MC: Jasmonate- and salicylate-mediated plant defense responses to insect herbivores, pathogens and parasitic plants. Pest Manag Sci. 2009, 65: 497-503. 10.1002/ps.1714.

    Article  CAS  PubMed  Google Scholar 

  46. Saito S, Hirai N, Matsumoto C, Ohigashi H, Ohta D, Sakata K, Mizutani M: Arabidopsis CYP707As Encode (1)-Abscisic Acid 8-Hydroxylase, a key enzyme in the oxidative catabolism of abscisic acid. Plant Physiol. 2004, 134: 1439-1449. 10.1104/pp.103.037614.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  47. An SH, Sohn KH, Choi HW, Hwang IS, Lee SC, Hwang BK: Pepper pectin methylesterase inhibitor protein CaPMEI1 is required for antifungal activity, basal disease resistance and abiotic stress tolerance. Planta. 2008, 228: 61-78. 10.1007/s00425-008-0719-z.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  48. Kogovsek P, Pompe-Novaka M, Baeblera S, Rottera A, Gowb L, Grudena K, Fosterb GD, Boonhamc N, Ravnikar M: Aggressive and mild Potato virus Y isolates trigger different specific responses in susceptible potato plants. Plant Pathol. 2010, 59: 1121-1132. 10.1111/j.1365-3059.2010.02340.x.

    Article  CAS  Google Scholar 

  49. Umate P: Genome-wide analysis of lipoxygenase gene family in Arabidopsis and rice. Plant Signal Behav. 2011, 3: 335-338.

    Article  Google Scholar 

  50. Xu L, Zhu L, Tu L, Liu L, Yuan D, Jin L, Long L, Zhang X: Lignin metabolism has a central role in the resistance of cotton to the wilt fungus Verticillium dahliae as revealed by RNA-Seq-dependent transcriptional analysis and histochemistry. J Exp Bot. 2011, 23: 1-15.

    CAS  Google Scholar 

  51. Santos RC, Marcellino LH, Monnerat RG, Gander ES: Mechanical damage in cotton buds caused by the boll weevil. Pesq Agropec Bras Brasília. 2003, 38: 1351-1356. 10.1590/S0100-204X2003001100015.

    Article  Google Scholar 

  52. Liebrand TWH, van den Berga GCM, Zhanga Z, Smita P, Cordewenerb JHG, America AHP, Sklenard J, Jonesd AME, Tamelinga WIL, Robatzekd S, Thommaa BPHJ, Joostena MHAJ: Receptor-like kinase SOBIR1/EVR interacts with receptor-like proteins in plant immunity against fungal infection. Proc Natl Acad Sci U S A. 2013, 24: 10010-10015.

    Article  Google Scholar 

  53. Liu J, Elmore LM, Lin ZD, Coaker G: A receptor-like cytoplasmic kinase phosphorylates the host target RIN4, leading to the activation of a plant innate immune receptor. Cell Host Microbe. 2011, 2: 137-146.

    Article  Google Scholar 

  54. Son GH, Wan J, Kim HJ, Nguyen XC, Chung WS, Hong JC, Stacey G: Ethylene-Responsive Element-Binding Factor 5, ERF5, Is involved in chitin-induced innate immunity response. Mol Plant Microbe Interact. 2012, 25: 48-60. 10.1094/MPMI-06-11-0165.

    Article  CAS  PubMed  Google Scholar 

  55. McGrath KC, Dombrecht B, Manners JM, Schenk PM, Edgar CI, Maclean DJ, Scheible W, Udvardi MK, Kazan K: Repressor- and activator-type ethylene response factors functioning in jasmonate signaling and disease resistance identified via a genome-wide screen of arabidopsis transcription factor gene expression. Plant Physiol. 2005, 139: 949-959. 10.1104/pp.105.068544.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  56. Hu Y, Donga Q, Yua D: Arabidopsis WRKY46 coordinates with WRKY70 and WRKY53 in basal resistance against pathogen Pseudomonas syringae. Plant Sci. 2012, 185: 288-297.

    Article  PubMed  Google Scholar 

  57. Rasmussen MW, Roux M, Petersen M, Mundy J: MAP kinase cascades in Arabidopsis innate immunity. Front Plant Sci. 2012, 24: 169-

    Google Scholar 

  58. Asai T, Tena G, Plotnikova J, Willmann MR, Chiu W, Gomez-Gomez L, Boller T, Ausubel FM, Sheen J: MAP kinase signalling cascade in Arabidopsis innate immunity. Nature. 2001, 415: 977-983.

    Article  Google Scholar 

  59. Howe GA, Jander G: Plant immunity to insect herbivores. Annu Rev Plant Biol. 2008, 59: 41-66. 10.1146/annurev.arplant.59.032607.092825.

    Article  CAS  PubMed  Google Scholar 

  60. Torres MA, Dangl JL, Jones JDG: Arabidopsis gp91phox homologues AtrbohD and AtrbohF are required for accumulation of reactive oxygen intermediates in the plant defense response. Proc Natl Acad Sci U S A. 2002, 99: 517-522. 10.1073/pnas.012452499.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  61. Leitner M, Boland W, Mithöfer A: Direct and indirect defences induced by piercing-sucking and chewing herbivores in Medicago truncatula. New Phytol. 2005, 167: 597-606. 10.1111/j.1469-8137.2005.01426.x.

    Article  CAS  PubMed  Google Scholar 

  62. Bostock RM: Signal crosstalk and induced resistance: straddling the line between cost and benefit. Annu Rev Phytopathol. 2005, 43: 545-580. 10.1146/annurev.phyto.41.052002.095505.

    Article  CAS  PubMed  Google Scholar 

  63. O'Donnell PJ, Calvert C, Atzorn R, Wasternack C, Leyser HM, Bowles DJ: Ethylene as a signal mediating the wound response of tomato plants. Science. 1996, 274: 1914-1917. 10.1126/science.274.5294.1914.

    Article  PubMed  Google Scholar 

  64. Laudert D, Weiler EW: Allene oxide synthase: a major control point in Arabidopsis thaliana octadecanoid signalling. Plant J. 1998, 5: 675-684.

    Article  Google Scholar 

  65. Schmelz EA, Alborn HT, Banchio E, Tumlinson JH: Quantitative relationships between induced jasmonic acid levels and volatile emission in Zea mays during Spodoptera exigua herbivory. Planta. 2003, 216: 665-673.

    CAS  PubMed  Google Scholar 

  66. von Dahl CC, Winz RA, Halitschke R, Kuhnemann F, Gase K, Baldwin IT: Tuning the herbivore-induced ethylene burst: the role of transcript accumulation and ethylene perception in Nicotiana attenuata. Plant J. 2007, 51: 293-307. 10.1111/j.1365-313X.2007.03142.x.

    Article  CAS  PubMed  Google Scholar 

  67. Liu Y, Schiff M, Serino G, Deng X-W, Dinesh-Kumar SP: Role of SCF Ubiquitin-Ligase and the COP9 Signalosome in the N Gene–Mediated Resistance Response to Tobacco mosaic virus. Plant Cell. 2002, 14: 1483-1496. 10.1105/tpc.002493.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  68. Yang C-W, Lamothe RG, Ewan RA, Rowland O, Yoshioka H, Shenton M, Ye H, O’Donnell E, Jones JDG, Sadanandoma A: The E3 ubiquitin ligase activity of arabidopsis PLANT U-BOX17 and its functional tobacco homolog ACRE276 are required for cell death and defense. Plant Cell. 2006, 18: 1084-1098. 10.1105/tpc.105.039198.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  69. Trujillo M, Ichimura K, Casais C, Shirasu K: Negative regulation of PAMP-triggered immunity by an E3 ubiquitin ligase triplet in Arabidopsis. Curr Biol. 2008, 18: 1396-1401. 10.1016/j.cub.2008.07.085.

    Article  CAS  PubMed  Google Scholar 

  70. Hubert DA, Tornero P, Belkhadir Y, Krishna P, Takahashi A, Shirasu K, Dangl JL: Cytosolic HSP90 associates with and modulates the Arabidopsis RPM1 disease resistance protein. EMBO J. 2003, 22: 5679-5689. 10.1093/emboj/cdg547.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  71. Takahashi A, Casais C, Ichimura K, Shirasu K: HSP90 interacts with RAR1 and SGT1 and is essential for RPS2-mediated disease resistance in Arabidopsis. Proc Natl Acad Sci U S A. 2003, 100: 11777-11782. 10.1073/pnas.2033934100.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  72. Seo Y-S, Lee S-K, Song M-Y, Suh J-P, Hahn T-R, Ronald P, Jeon J-S: The HSP90-SGT1-RAR1 molecular chaperone complex: a core modulator in plant immunity. J Plant Biol. 2008, 51: 1-10. 10.1007/BF03030734.

    Article  CAS  Google Scholar 

  73. Lu R, Malcuit I, Moffett P, Ruiz MT, Peart J, Wu AJ, Rathjen JP, Bendahmane A, Day L, Baulcombe DC: High throughput virus-induced gene silencing implicates heat shock protein 90 in plant disease resistance. EMBO J. 2003, 22: 5690-5699. 10.1093/emboj/cdg546.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  74. Liu Y, Burch-Smith T, Schiff M, Feng S, Dinesh-Kumar SP: Molecular chaperone Hsp90 associates with resistance protein N and its signaling proteins SGT1 and Rar1 to modulate an innate immune response in plants. J Biol Chem. 2004, 279: 2101-2108. 10.1074/jbc.M310029200.

    Article  CAS  PubMed  Google Scholar 

  75. Monnerat RG, Dias SC, Oliveira-Neto OB, Nobre SD, Silva-Werneck JO, Grossi de Sá MF: Comunicado Técnico. Embrapa Recursos Genéticos e Biotecnologia. Criação massal do bicudo do algodoeiro Anthonomus grandis em laboratório. 2000

    Google Scholar 

  76. Schroeder A, Mueller O, Stocker S, Salowsky R, Leiber M, Gassmann M, Lightfoot S, Menzel W, Granzow M, Ragg T: The RIN: an RNA integrity number for assigning integrity values to RNA measurements. BMC Mol Biol. 2006, 7: 3-10.1186/1471-2199-7-3.

    Article  PubMed Central  PubMed  Google Scholar 

  77. Xie F, Sun G, Stiller JW, Zhang B: Genome-wide functional analysis of the cotton transcriptome by creating an integrated EST database. Plos One. 2011, 6: 1-12.

    Article  Google Scholar 

  78. Li H, Durbin R: Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009, 25: 1754-1760. 10.1093/bioinformatics/btp324.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  79. Wang L, Feng Z, Wang X, Wang X, Zhang X: DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010, 26: 136-138. 10.1093/bioinformatics/btp612.

    Article  PubMed  Google Scholar 

  80. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M: Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2001, 21: 3674-3676.

    Article  Google Scholar 

  81. Alexa A, Rahnenfuhrer J, Lengauer T: Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006, 22: 1600-1607. 10.1093/bioinformatics/btl140.

    Article  CAS  PubMed  Google Scholar 

  82. Thimm O, Bläsing O, Gibon Y, Nagel A, Meyer S, Krüger P, Selbig J, Müller LA, Rhee SY, Stitt M: MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J. 2004, 6: 914-939.

    Article  Google Scholar 

  83. Rozen S, Skaletsky H: Primer3 on the WWW for general users and for biologist programmers. Methods Mol Biol. 2000, 132: 365-386.

    CAS  PubMed  Google Scholar 

  84. Scanlon MJ, Ohtsu K, Timmermans MC, Schnable PS: Laser microdissection-mediated isolation and in vitro transcriptional amplification of plant RNA. Curr Protoc Mol Biol. 2009, Chapter 25: Unit 25A.3-

    PubMed  Google Scholar 

  85. Zhao S, Fernald RD: Comprehensive algorithm for quantitative real-time polymerase chain reaction. J Comput Biol. 2005, 12: 1047-1064. 10.1089/cmb.2005.12.1047.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  86. Artico S, Nardeli SM, Brilhante O, Grossi-de-Sa MF, Alves-Ferreira M: Identification and evaluation of new reference genes in Gossypium hirsutum for accurate normalization of real-time quantitative RT-PCR data. BMC Plant Biol. 2010, 10: 49-10.1186/1471-2229-10-49.

    Article  PubMed Central  PubMed  Google Scholar 

  87. Pfaffl MW, Horgan GW, Dempfle L: Relative expression software tool (REST©) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002, 30: 2-10. 10.1093/nar/30.2.e2.

    Article  Google Scholar 

  88. Bailey TL, Williams N, Misleh C, Li WW: MEME: discovering and analyzing DNA and protein sequence motifs. Nucleic Acids Res. 2006, 34: W369-W373. 10.1093/nar/gkl198.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  89. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S: MEGA5: Molecular Evolutionary Genetics Analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011, 10: 2731-2739.

    Article  Google Scholar 

  90. Whelan S, Goldman N: A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach. Mol Biol Evol. 2001, 5: 691-699.

    Article  Google Scholar 

  91. Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O: New algorithms and methods to estimate maximum-likelihood phylogenies: as sessing the performance of PhyML 3.0. Syst Biol. 2010, 59 (3): 307-321. 10.1093/sysbio/syq010.

    Article  CAS  PubMed  Google Scholar 

Download references


We thank Prof. S.M. Tsai, Cell and Molecular Biology Laboratory, CENA/USP, for maintaining the multiuser LMD facility (FAPESP 2009/53998-3). APM acknowledges CNPq (Grant 310.612/2011-0). The authors gratefully thank Marco A. Takita and Carolina M. Rodrigues for use of the computational structure in the Laboratory of Biotechnology, Center for Citrus Sylvio Moreira, Cordeirópolis-Sao Paulo, Brazil. We thank Julia Lambret-Frotté for help us with the phylogenetic analysis. This work was part of S.A’s Ph. D. research in Genetics at the Department of Genetics of the Universidade Federal do Rio de Janeiro (UFRJ) and was supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq; M. Alves-Ferreira: # 306025/2010-8) and Fundação de Amparo à Pesquisa do Rio de Janeiro (FAPERJ).

Author information

Authors and Affiliations


Corresponding author

Correspondence to Marcio Alves-Ferreira.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

MA-F, MFG-S and SA planned and supervised the study. MA-F and SA contributed to the design, execution of the experiments and drafted the manuscript. OBON and LLPM conducted inoculation and collection of samples. MRA and SA contributed to the RNA-seq analysis, categorization and annotations of differentially expressed transcripts. APM and SS contributed with the microdissection experiment. MA-F, SA contributed to the interpretation of the data and provided intellectual input. APM, OBON revised the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material


Additional file 1: JPEG file showing contig size distribution from the aligned reads of transcriptome sequencing of cotton flower buds.(TIFF 11 MB)


Additional file 2: Full list of the up-regulated and down-regulated cotton genes found in the RNA-seq analysis. Only those differentially expressed genes with significant expression changes (adjusted p-value ≤ 0.05, |logFC| ≥ 2.0) are shown. Columns A-G are the results from Blastx with Arabidopsis, and columns H-P are the results from the Blast2GO program. (XLS 780 KB)


Additional file 3: Molecular function (A) and Cellular component (B) ontologies overrepresented by the Gene Set Enrichment Analysis (GSEA).(TIFF 521 KB)


Additional file 4: Subset of differentially expressed genes (DEG) in cotton flower buds in response to cotton boll weevil larvae feeding that were annotated with the GO term “plasma membrane”.(XLS 44 KB)


Additional file 5: Subset of differentially expressed genes (DEG) in cotton flower buds in response to cotton boll weevil larvae feeding that were annotated with the GO term “cell wall”.(XLS 38 KB)

Additional file 6: Biological processes (BP) overrepresented by the Gene Set Enrichment Analysis (GSEA).(XLS 58 KB)


Additional file 7: The twenty-one cotton WRKY genes that were differentially expressed in our RNAseq experiment, the domain sequence used in the phylogenetic tree, the protein sequence and the nucleotide sequence.(XLS 77 KB)


Additional file 8: Comparison of qPCR and RNA sequencing expression data. Thirty-two differentially expressed genes responding to cotton boll weevil larvae feeding were analysed for RNA abundance using qPCR. The log fold-change in RNA sequencing was estimated and tested by DEGseq software, in which the count data were normalised to the total number of counts, taking the variance and the mean of the biological replicates into account for each library with reference to the control (blue bars). qPCR results (red bars) were calculated with the relative REST software using a mathematical model based on the PCR efficiencies and the mean CP deviation between the sample and control group of target genes and were normalised to the mean CP deviation of reference genes. We calculated the Pearson correlation coefficient (r) between the different methods for all transcripts. The correlation coefficient was 0.9406. (TIFF 8 MB)


Additional file 9: Primers used for the qPCR expression analysis of DEG in tissues isolated from cotton flower buds infested by cotton boll weevil larvae.(XLS 52 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit

The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Artico, S., Ribeiro-Alves, M., Oliveira-Neto, O.B. et al. Transcriptome analysis of Gossypium hirsutum flower buds infested by cotton boll weevil (Anthonomus grandis) larvae. BMC Genomics 15, 854 (2014).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: