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

Background 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. Results 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). Conclusion 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. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-854) contains supplementary material, which is available to authorized users.


Background
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,[14][15][16][17][18]. Highthroughput 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 tissuechewing 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).
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][40][41][42][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 downregulated 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 Blas-t2GO 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 pathogenderived 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).
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-, herbivoreand 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).

Main processes or pathways affected by the response to infection with cotton boll weevil larvae
Several subsets of genes were identified, including protein kinases, Ca 2+ -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, Ca 2+ -binding proteins and protein degradation An appropriate defence response to a biotic threat initially requires recognition of the threat. Herbivoreassociated 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" (  Only DEG with the GO term for the biological process "response to chitin" (GO:0010200) and/or "death" (GO:0016265) are shown. The complete list of DEG is shown in Additional file 2. *Genes present in both biological processes: response to chitin and death. Genes in bold were tested by qPCR. 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 Ca 2+ dependent manner. Mechanical wounding and insect attack trigger a transient increase of cytosolic Ca 2+ levels within minutes. This Ca 2+ oscillation acts as a mediator for stimulus responses to several Ca 2+ -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 Ca 2+ influx into the cytosol (Table 3) [30]. Targeted protein degradation via the ubiquitin/26Sproteasome 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 helixloop-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 [32][33][34].
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).

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 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.
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 upregulated 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).
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 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 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 C P values by a Pair-Wise Fixed Reallocation Randomization Test was performed. 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 germinlike protein (GER3) with oxalate oxidase activity leading to hydrogen peroxide (H 2 O 2 ) production, which is a ROS that is implicated in the herbivory-induced response in plants [17].

Discussion
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 leucinerich repeats and an intracellular kinase domain for the activation of downstream signalling, as well as receptorlike 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. GhWRKY33like, 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 (Ca 2+ ) 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 Ca 2+ levels [59]. For instance, the Egyptian cotton worm (Spodoptera littoralis) caused a transient increase in cytosolic Ca 2+ in Phaseolus lunatus cells adjacent to the insect damage [59]. Cotton genes encoding calmodulinlike 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 Ca 2+ 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 (O 2− ), hydrogen peroxide (H 2 O 2 ), singlet oxygen ( 1 O 2 ), 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 ubiquitinligase 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 PAMPtriggered 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 (H 2 O 2 ) 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 [70][71][72]. 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 downregulated, suggesting that larvae attempt to overcome the plant immune defence triggered by the effector.

Conclusion
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); Ca 2+ 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, herbivoryinduced 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. grandis infestation 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). Threemonth 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, basecalling, 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; http://www.leonxie.com) [77] using the BWA package [78]. Differential expression was estimated and tested with the DEGseq software package in Rbioconductor 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 log 2 (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 (http://www.blast2go.com/b2ghome) [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].

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 OD 260 . 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 C P deviation between the sample and control group of target genes and is normalised to the mean C P deviation of the reference genes [83]. The advantage of REST is the use of a subsequent statistical test of the analysed C P values by a Pair-Wise Fixed Reallocation Randomization Test [87].