The transcriptional response to the olive fruit fly (Bactrocera oleae) reveals extended differences between tolerant and susceptible olive (Olea europaea L.) varieties

The olive fruit fly Bactrocera oleae (Diptera: Tephritidae) is the most devastating pest of cultivated olive (Olea europaea L.). Intraspecific variation in plant resistance to B. oleae has been described only at phenotypic level. In this work, we used a transcriptomic approach to study the molecular response to the olive fruit fly in two olive cultivars with contrasting level of susceptibility. Using next-generation pyrosequencing, we first generated a catalogue of more than 80,000 sequences expressed in drupes from approximately 700k reads. The assembled sequences were used to develop a microarray layout with over 60,000 olive-specific probes. The differential gene expression analysis between infested (i.e. with II or III instar larvae) and control drupes indicated a significant intraspecific variation between the more tolerant and susceptible cultivar. Around 2500 genes were differentially regulated in infested drupes of the tolerant variety. The GO annotation of the differentially expressed genes implies that the inducible resistance to the olive fruit fly involves a number of biological functions, cellular processes and metabolic pathways, including those with a known role in defence, oxidative stress responses, cellular structure, hormone signalling, and primary and secondary metabolism. The difference in the induced transcriptional changes between the cultivars suggests a strong genetic role in the olive inducible defence, which can ultimately lead to the discovery of factors associated with a higher level of tolerance to B. oleae.


Introduction
The olive fruit fly, Bactrocera oleae (Diptera: Tephritidae), is arguably the single largest threat of cultivated olive (Olea europaea L.). Since olive domestication, B. oleae is the most devastating pest in the Mediterranean basin and the recent invasion of California suggests that the fruit fly follows its host [1,2]. Phytophagous larvae feed exclusively on olive fruits [3]. Depending on climate conditions, the fruit fly can produce severe crop damage and significant economic loss in the whole olive sector [4].
Adult females pierce the olive and lay eggs under the exocarp. Hatched larvae progressively consume the olive pulp, causing tissue loss, premature fruit drop and reduction of oil yield. Moreover, fly infestation increases olive oil acidity and peroxide value, as well as musty and earthy off-flavours, extensively reducing oil quality (e.g. downgrading extra virgin olive oil to less valuable categories). Indirect effects are mainly due the presence of necrotic areas and microorganisms in feeding tunnels [5][6][7][8][9]. The conventional B. oleae management relies on chemical insecticides and/or traps [4]. The olive fly, as many other pests, can acquire resistance to pesticides [10], increasing the need for more effective biological or integrated control methods [11][12][13].
Considering the unrivalled olive oil-health benefits, research on the olive-fruit fly interaction has the long-term potential to influence not only the olive oil production but also the health-promoting properties of the olive oil. The economic impact of the B. oleae does not match the research efforts aimed to examine the interaction between the olive and its key enemy. For instance, compared to other biotic stresses, very few studies shed light on the mechanisms underlying olive defence and resistance at molecular level [14,15].
Olive varieties can have different levels of tolerance to the olive fruit fly [4,16]. In a comparison of different cultivars, the percentage of infestation ranged from less than 10% to up to 31% [16]. This difference is significant because a percentage of exit holes lower than 10% is compatible with a production of high quality olive oil in absence of chemical control methods, if adequate harvesting and storing procedures are followed [7]. A different tolerance to the olive fruit fly is evident not only comparing cultivars of different origins but also analysing regional germplasm [16,17]. For instance, among the twenty varieties that constitute the olive germplasm of the Campania region (Southern Italy) [18], the 'Ortice' and 'Ruveia' are reported to be highly susceptible and tolerant, respectively [19].
The basis of the different tolerance to the olive fruit fly is expected to be complex [20]. It may rely on mechanical barriers (e.g. aliphatic waxes), chemical factors (e.g. oleuropein, cyanidine), morphological characteristics (e.g fruit size) and their combination. Similarly, the relative prominence and the contribution of these factors are yet to be fully clarified [21][22][23][24][25]. Moreover, the mechanism underlying the different levels of tolerance to the B. oleae has never been studied at molecular level. Finally, it is not known whether those features are constitutively expressed or induced by the fruit fly feeding.
To address these points, we studied the molecular response of the drupe in two olive cultivars with different levels of tolerance to the fruit fly. To this aim, we first generated a catalogue of more than 80.000 unigenes by next-generation pyrosequencing and then developed a microarray layout as affordable tool for functional genomics in olive. Our study provides insights into the molecular reaction of the drupe to larva feeding and illustrates the complexity and the differences of the drupe defence response of varieties with different levels of fruit fly susceptibility.

Plant material
Two olive varieties were studied because of their contrasting tolerance to B. oleae, 'Ortice' (susceptible to the olive fruit fly) and 'Ruveia' (tolerant) [19]. Olives were harvested from field plants at the Azienda Agricola Regionale Sperimentale "Improsta" (Eboli, Salerno) when at the Jaen Ripening Index (JRI) 2 [26]. The percentage of olive fly attack was calculated on 300 drupes per variety (10 replicated groups of 30 randomly collected drupes). Statistical differences were evaluated with a Student t-test. For molecular studies, olives with visible symptoms of pathogen attack as well as those with fly exit holes were discarded. Olives at different attack stages (punctured, with II or III instar larva, or pupae) were pooled and used with undamaged olives for massive parallel sequencing, in order to create a comprehensive repertoire of genes modulated by B. oleae. For the study of differential gene expression by microarray analysis, we used for each cultivar drupes with feeding II or III instar larvae and undamaged controls. For massive sequencing and for microarray analysis, olives were sliced under a light microscope to remove the larva. At this stage, drupes with sterile stings and dead/inactive larvae were also discarded. Slices were frozen in liquid nitrogen and stored at -80˚C until use.

sequencing
Total RNA was isolated from 200 mg of drupe with the RNeasy Plant Mini Kit (Qiagen Hilden, Germany) and treated with DNase I (Ambion, Austin, TX, USA). Three replicates, each consisting of five olives, were used for each biological condition (punctured olives, olives with feeding II or III instar larva, olives with pupae and undamaged olives). cDNAs were synthesized using the SMART PCR cDNA Synthesis kit (Clontech, Palo Alto, CA, USA). First strand synthesis was performed using eight μg of total RNA as described [27]. Double stranded cDNAs were purified using the QIAquick PCR purification kit (Qiagen, Hilden, Germany) and quantified with a fluorimeter (Victor 2, Perkin Elmer, Wellesley, MA, USA). To estimate cDNA quality and fragment length, samples were separated on a 1.5% agarose gel. cDNA from olives characterized by different stages of fly attack were pooled together to obtain a representative 'damaged' sample. Two cDNA libraries (damaged and undamaged olives) per cultivar ('Ortice' and 'Ruveia') were sequenced. Approximately five μg of cDNA were sheared into small fragments via nebulization. The four shotgun cDNA libraries were sequenced using a 454 GS FLX+ Titanium Sequencer (Roche Diagnostics Corporation, Basel, Switzerland). The adaptor-trimmed 454 reads were assembled using the GSAssembler Software (Roche Diagnostics Corporation, Basel, Switzerland). To annotate unigenes (contigs and singletons), a blastXbased similarity search (e-value 1e-5) was used, querying the NCBI non-redundant (nr) database. We mapped the GI identifiers (http://www.ncbi.nlm.nih.gov/) of the best blastX hits to the UniprotKB protein database (http://www.uniprot.org/) in order to extract Gene Ontology (GO, http://www.geneontology.org/) and KEGG orthology (KO, http://www.genome.jp/ kegg/) terms.

Microarray preparation and hybridization
CustomArray™ (CustomArray Inc., Bothell, WA, USA) microarrays containing over 90k probes were built at the ENEA-Trisaia Research Center (Rotondella, MT, Italy). Probes were designed using the ProbeWeaver software (CustomArray Inc., Bothell, WA, USA), based on the 454 pyrosequencing results of four cDNA libraries obtained from damaged and undamaged drupes of the 'Ortice' and 'Ruveia' varieties.
The chip layout consists of 61,825 olive probes (60,706 non-redundant) out of the 87,720 sequences of the pooled library. The layout includes around 20.000 quality control spots. Probes were 35-to 40-mers with a melting temperature of 70-75˚C and were synthesized on microarrays through the CustomArray Synthesizer™ (CustomArray Inc., Bothell, WA, USA).
Total RNA was isolated from 200 mg of finely ground olive tissue. A purification step was performed in a 1:1:2 (v/v/v) solution of phenol, chloroform and RNA Extraction Buffer (1 M Tris-HCl pH 8.5, 5 M NaCl, 0.5 M EDTA pH 8.0, 10% SDS) twice, followed by a chloroform purification. RNA isolation from the aqueous phase and quality control were performed as described [28]. Samples with a 260/280 nm absorbance ratio higher than 1.8 and a 260/230 nm absorbance ratio higher than 2.0 were used for subsequent experiments. For each biological condition, three independent samples (i.e. from different trees) were obtained as pools of three to five independently extracted technical replicates. Two μg of total RNA were retrotranscribed using the RNA ampULSe Amplification and Labeling Kit for CustomArray™ microarrays (Kreatech Biotechnology). cRNA labeling was carried out with the Cy5-ULS (Cyanine-Universal Linkage System). Unincorporated dye was removed with a KREApure purification column (Agilent). Labeling yield and quality were assessed using a NanoDrop 1000 (Thermo Scientific). cRNAs were fragmented for 20 min at 95˚C in 200 mM Tris-Acetate (pH 8.1), 500 mM potassium acetate and 150 mM molybdenum acetate. The fragmented Cy5-cRNAs were added to the hybridization solution (6X SSPE, 0.05% Tween-20, 20mM EDTA, 25% formamide, 100 ng/μl salmon sperm DNA, 0.04% SDS) and poured in the hybridization chambers of the prehybridized microarrays, following the CustomArray™ Hybridization protocol. After an overnight incubation at 45˚C in a rotating rotisserie oven, microarrays were treated with six washing steps: firstly, in 6X SSPE and 0.05% Tween-20 for 5' at 45˚C, then with 3X and 0.5X SSPE and 0.05% Tween-20 for 1' at room temperature, with 2X PBS and 0.1% Tween-20, and finally twice with 2X PBS for 1' at room temperature.
For the validation of the microarray results, Real time PCR was carried out as described [14]. Primers and their main features are reported in S1 Table. Microarray data analyses Microarray slides were scanned with the GenePix Pro microarray scanner (Axon Instruments) and data processed with the CustomArray™ Microarray Imager software (CustomArray Inc., Bothell, WA, USA). A maximum threshold of 0.20 for the coefficient of variance (of the spots corresponding to identical probes) was applied to control intra-chip hybridization variability. The correlation among the three technical replicates for each experimental condition was assured by a minimum Pearson coefficient (R) of 0.99. Stripping and preparation of the microarrays for re-hybridization was performed twice, considering the three technical replicates, according to the manufacturer's instruction (CustomArray Inc., Bothell, WA, USA). Raw values were normalized on the median of the intensities by the ProbeWeaver software (Custo-mArray Inc., Bothell, WA, USA) using the quantile algorithm. Pairwise analysis of differential expression was assessed with a Welch's t-test, followed by a Bonferroni-Hochberg FDR adjustment with a cut off p-value of 0.05, using R [29]. To filter out weakly expressed sequences, we calculated the average and standard deviation of the expression value of all empty and negative controls and set as threshold the mean value plus two times the standard deviation [30]. Of the filtered, significantly differentially expressed probes, only those with greater than 2-fold increase or 2-fold decrease in expression compared to the control condition were used for further analysis. Principal Component Analysis was performed in R [29]. K-means clustering analysis of the differentially expressed genes (DEGs) was carried out with the Multi Experiment Viewer (MeV) software [31]. Dataset grouping was performed with the average linkage clustering methods based on Euclidean distances. Different number of clusters (k fixed from 8 to 16) were compared and k = 10 was selected to maximize partitioning and to avoid empty or very small groups. The corresponding expression graphs were visualized with MeV. We followed the Minimum Information About a Microarray Experiment (MIAME) guidelines for microarray analysis and verification [32].

Functional annotation of differentially expressed genes
Functional annotation was carried out by sequence analysis using the Blast2GO software [33]. Briefly, for each sequence corresponding to a differentially expressed probe, a blastX similarity search (e-value < 1e-6) against the non-redundant (nr) NCBI database (Non-redundant Gen-Bank CDS translations including RefSeq, PDB, SwissProt, IR and PRF) was performed to retrieve a maximum of 20 top homologous hits per query. The GO-term mapping and annotation were retrieved using NCBI as well as non-redundant reference protein database (PSD, UniProt, RefSeq, GenPept, PDB Full Gene Ontology DB). Sequences with a blast hit that could not be mapped and annotated were then blasted (blastX) against the Arabidopsis thaliana protein sequences and the Oryza sativa protein sequences database. Additional annotations (e.g. the recovery of implicit "Biological Process" and "Cellular Component" GO-terms from "Molecular Function" annotations) were implemented using ANNEX 5.0 [34]. Completion of the functional annotation with protein domain information was performed with InterProScan 5.0. A plant GO-Slim reduction was carried out to summarize the functional content of the dataset. To complete the functional annotation Sma3s was also used [35]. Plant taxonomic division from UniProt was selected as source of annotations, and the transcriptomic sequences were enriched with keywords, InterPro domains and motifs, and GO terms. We generated multi-pie chart summary to present meaningful GO identifiers that are at hand yet not excessively general. When possible, sequences that could not be mapped were named (but not annotated) after a tblastX search against the nr NCBI database to identify sequences similar to the query based on their coding potential.

Results
The already reported different susceptibility to B. oleae of the two varieties was evaluated by measuring the percentage of attack at the time of harvesting. The susceptible variety 'Ortice' had an attack index 2.5 times higher than the 'Ruveia' (p<0.01). At sample harvest, the infestation level was 13.3% in 'Ortice' and 5.3% in 'Ruveia'.
Pyrosequencing of the four cDNA libraries from damaged (pool of different stages of fly attack) and control undamaged fruits of 'Ortice' and 'Ruveia' generated 695,511 reads, with an average length of 361 bp (N50: 421 bp; N75: 359 bp) ( Table 1).
We assembled the raw reads from the four libraries together. More than 80% of raw sequences were included in the assembly and 72,662 remained as singleton ( Table 2).
The assembly yielded 15,058 contigs, with an average length of 884 bp and 78% of the sequences longer than 500 bp (Panel A in S1 Fig). Approximately 40% of the contigs were composed of 2-10 reads. The majority of the contigs (28%) were low-coverage (6-10 reads) and around 5% were high-coverage (> 100 reads). The unigenes list and their annotation is reported in S2 Table. An overview of the putative functions of the unigenes, based on a blastXsimilarity search annotation, is presented in S2 To identify main biological pathways, we mapped the annotations to reference canonical pathways in the KEGG database by using KO identifiers (i.e. a classification of orthologous genes defined by KEGG). The 5,355 KO allocated genes were assigned to 265 KEGG categories. The most abundant processes/metabolic pathways were "biosynthesis of secondary metabolites" (185 unigenes); "microbial metabolism in diverse environments" (89); "spliceosome" (74); "biosynthesis of plant hormones" (72); "ribosome" (61); "RNA transport" (61) and "protein processing in endoplasmic reticulum" (60)

Differential gene expression
To profile the variation in gene expression of olives during B. oleae feeding and to characterize the differences between the two varieties with contrasting susceptibility to the olive fruit fly, unigenes were used to build a CustomArray™ (CombiMatrix Corporation). The chip layout consists of 61,825 olive probes out of the 87,720 assembled sequences (Panel B in S1 Fig). The DNA microarrays were then hybridized with labeled cRNA prepared from undamaged olives of the varieties 'Ortice' or 'Ruveia' and corresponding samples with feeding larvae.
After normalization (S3 Fig), the expression data relative to all the olive probes were analysed by PCA to explore the variation among the different conditions based on gene expression states (Fig 1).
The variability among the different conditions and biological replicates could be summarized and visualized in two dimensions. The first component, which accounted for 46.0% of total variance, appears to represent the overall induction due to the fruit fly, as it discriminates the different biological response of the two varieties. The second component, which explained 14.1% of the variance, mainly discriminated the biological replicates for each condition. Because of a clear B. oleae effect, the limited difference among the control conditions of the two cultivars and the dispersion of the biological replicates, we retained all conditions and replicates for further analysis.
To identify differentially expressed probes, we applied the following selection criteria: a pvalue <0.05 (Welch t-test, followed by a Bonferroni-Holchberg FDR correction for multiple testing) and a |log 2 Ratio|> 1. The number of differentially expressed genes, after removal of duplicates, is reported in Table 3.
DEGs and their annotation are listed in S1 File. The heatmap illustrates a weak linkage among the four conditions (S4 Fig). To validate the microarray results, the expression of six genes (four differentially expressed genes and two non-affected by B. oleae) was analyzed by Real Time-PCR. The results were consistent to the microarray data (S5 Fig).  The smallest difference, in terms of differentially expressed genes, was present between the undamaged olives of the two cultivars. The response to the B. oleae feeding of the susceptible cultivar 'Ortice' also involved a limited set of genes and the majority (77%) were down regulated. The response of the more tolerant variety 'Ruveia' involved the highest number of both up-regulated and down-regulated genes. When attacked, olives of the 'Ruveia' cultivar differentially expressed more than 1000 genes compared to attacked drupes of the 'Ortice'. To investigate the commonality between the olive varieties and their response, DEGs in the pairwise comparisons were matched. The response of 'Ruveia' is highly specific, as indicated by the limited number of common genes. Around half of the genes down-regulated in 'Ortice' were also down-regulated in 'Ruveia' (Fig 2).
The majority of the genes overexpressed in 'Ruveia' after B. oleae attack were also upregulated when comparing the damaged 'Ruveia' and 'Ortice' drupes. Overall, the comparison of the various response to the B. oleae indicated that two varieties have a markedly different reaction.

Transcript clustering
Considering the limited overlap of the molecular response of the two varieties, the expression level of the statistically significant genes in at least one of the four comparisons was processed by K-means clustering analysis to infer possible co-expressed or co-regulated genes (Fig 3).
This analysis allowed to group the 2848 DEGs in 10 clusters (S3 Table). K-means clustering indicated that it was not possible to identify a cluster in which genes are overexpressed in the two attacked varieties with a similar pattern, further indicating a very limited overlap between the responses of the two olive varieties. The GO annotation (level 4) in the "Biological Process" domain of the genes belonging to the clusters (sequence cut-off 5%) is reported in Fig 4. To ease the comparison, the graph reports, for each cluster, the different GO categories in relative terms. The annotation indicated that K7 and K4 are the most complex clusters because they include, respectively, 12 and 11 out of the 15 considered GO terms. Taking as reference the average expression pattern, clusters K7 and K4 mainly comprise genes that are, respectively, highly and mildly up-regulated exclusively in the damaged drupes of the more tolerant cultivar 'Ruveia'. Similarly, genes belonging to clusters K1 and K6, specifically modulated in 'Ruveia' drupes, were annotated with a high number of GO categories.
The most present GO terms in the different clusters were related to chemical reactions and pathways involving macromolecule, such as "macromolecule metabolic process" (present in all the ten clusters), "cellular macromolecule metabolic process" and "cellular nitrogen compound metabolic process" (present in nine clusters). K4 is the only cluster containing transcripts related to "carbohydrate metabolic processes" while the "post-embryonic development" term was exclusively present in K7. This analysis indicated the degree of specificity of the molecular events and biological processes characterizing the response of the tolerant cultivar 'Ruveia' to B. oleae attack.

Expression analysis of the 'Ruveia' response to fruit fly
We used GO terms association (obtained by a sequence similarity search) to elucidate the biological objective to which the DEGs contribute in the response to the B. oleae feeding. GO analysis indicated that the transcriptional reconfiguration involved a range of biological processes (S6 Fig). To provide an optimal view of the dataset's most relevant terms, a summary of the annotation results is presented as a multi-level pie that shows the lowest GO terms per branch that fulfil our annotation criterion weight (Fig 5).
The functional profile of both down and upregulated sequences indicated that the most abundant GO term categories were related to chemical reactions and pathways resulting in the formation of substances, including protein modification processes. Collectively, response to stress and to abiotic or biotic stimulus ranked the first most affected biological process for the down-regulated genes and at the eighth place for the up-regulated sequences.
Among the up-regulated genes involved in the response to stress, there were transcripts coding for proteins involved in signaling, including receptors and transcription factors, such as two leucine-rich repeat receptors, nine WD-repeat containing proteins, five WRKY transcription factors, and two rpm1-like proteins. The latter are essential regulator of plant defense Drupe response to the olive fruit fly and are typically associated to the resistance to pathogens [36]. Other up-regulated genes typically associated to pathogenesis included two PR proteins, one disease resistance response protein, one late blight resistance protein homologue and a NBS domain resistance protein.
Among the genes that can directly deter insect growth, there were three beta-glucosidases. In fruits, beta-glucosidase activity plays an important role in the oleuropein metabolism, catalyzing its hydrolysis into a toxic glutaraldehyde-like structure that acts as defense mechanism against insects [37]. Moreover, four serine carboxypeptidase-like (SCPL) genes were up-regulated. In tomato, Arabidopsis and rice, some SCPL proteins are wound inducible and connected to jasmonic acid (JA) pathway [38,39]. Other up-regulated genes possibly related to JA pathway include two phospholipases. Genes typically associated to abiotic stress were also upregulated, such as two transcripts coding for dehydration-induced proteins, two heat-shock proteins and a bobber1-like protein. The 'Ruveia' response to B. oleae also included genes involved in phytohormone signaling such as those coding for proteins involved in auxin (i.e. auxin binding protein, auxin efflux carrier family protein, auxin response factor, indole-3 acetic acid amino-acid hydrolase), brassinosteroid (i.e. brassinosteroid insensitive 1-associated receptor kinase 1 like, bri1 suppressor 1 like) and ethylene signaling (i.e. four ethylene responsive transcription factors). Besides stress hormones, the plant perception and the signal transduction of B. oleae infestation seems to involve also other stress-related cellular messengers such ROS and calcium, as indicated by the up-regulation of five calcium-dependent CBLinteracting serine/threonine protein kinases, a calcium transporting ATPase plasma membrane, four glutaredoxin, and one oxoglutarate genes. KEGG pathway map analysis of annotated enzymatic activities indicated that purine metabolism ranked first as number of sequences. Among other physiological processes, variation in purine metabolism associates to fruit ripening and to uracil salvage activity during senescence. Moreover, plants produce toxic secondary metabolites from pyrimidines that act as defense compounds [40]. The "Biosynthesis of antibiotics" KEGG reference pathway ranked first considering the number of enzymatic activities. This molecular interaction network diagram comprises an ample number of reactions that in plants are associated to terpenoid backbone biosynthesis, the shikimate pathway and the biosynthesis of secondary metabolites. Starch and sucrose metabolism ranked second for both the number of sequences and enzymatic activities.  To understand the inducible factors that can account for the different tolerance of the two varieties, we compared the gene expression levels between 'Ortice' and 'Ruveia' drupes when infested by B. oleae. The total number of DEGs was less than half compared with the response of the tolerant cultivar 'Ruveia' to the fruit fly. The proportion of the up-regulated genes (57%) was higher compared to the 'Ruveia' response, indicating that the different tolerance between the cultivars lies in the activation of a wide number of genes. In the Biological Process domain, differentially expressed sequences were assigned to a number of processes. At the GO level 2, the largest difference in ranking between the up-regulated and down-regulated sequences was for "multicellular organismal process". On the opposite, "metabolic process", in comparative terms, it is the most characterizing GO term in the down-regulated sequences (S7 Fig). A summary of the GO annotation results for the down-regulated or the up-regulated sequences is presented as a multi-level pie (Fig 6).
The dominant GO term for both upregulated and downregulated sequences was "singleorganism process". "Response to stimulus" ranked fourth for the up-regulated sequences and ninth for the down-regulated sequences.
Among the up-regulated genes, eight sequences were associated to the GO term "response to biotic stimulus" such as a WRKY transcription factor 17, two putative serine/threonine kinases, and a Serine carboxypeptidase like. Genes coding for proteins putatively involved in phytohormone signaling included an auxin efflux carrier family protein and abri1 suppressor 1 like. These genes were also upregulated in 'Ruveia' cultivar in response to the olive fruit fly. In total, almost 95% of the genes up-regulated in the comparison between the attacked drupes of the two cultivars were also upregulated in the damaged fruits of the tolerant variety ('Ruveia'). The data indicated that the different tolerance to the fruit fly in the two cultivars under investigation is mainly inducible.

Discussion
The use of olive varieties that are highly tolerant to the fruit fly is an important element to reduce economic loss and use of chemical pesticides. Knowledge on the molecular aspects underlying the different resistance response to the fruit fly can assist the screening of more suitable genotypes and ultimately, contribute to the development of new integrated control strategies. For these reasons, we compared two olive cultivars with contrasting resistance levels to B. oleae. To this aim, we first generated a large collection of unigenes from our NGS data. The assembly did not include about 20% of raw sequences, a limited number considering the abundance of short repeats and the high level of heterozygosity of the olive [27,41]. The reduced number of unassembled reads should have been also affected by the smaller number of cultivars and plant tissues under investigation in comparison with other works [41]. The 454 DNA sequencing method confirmed to be an effective technology in revealing the expression of a large number of genes in a non-model organism mainly because its longer read length [42].
We pooled the pyrosequencing samples derived from a number of conditions to generate a microarray based on the drupe sequence information and obtained a number of unigenes comparable to other works in olive [41]. Similarly, the percentage of annotation of our assembled dataset was comparable or higher to previous olive ESTs sequencing efforts [27,41,43]. The setup of a microarray allows analyses of gene expression and functional genomics studies in olive at a fraction of the cost of next generation sequencing. Moreover, microarray data are more computationally tractable and do not require an extended bioinformatics effort [44]. Our layout specifically focused on the drupe and was based on the widely used Combimatrix technology, hence extending currently available options [45,46].
A key interest in this study was to compare the response of genotypes with different level of tolerance to the olive fruit fly. In absence of infestation, very little difference was noticed between the two cultivars. This is justified considering that we analysed plants present in the same environment, not only to reduce environmental variability but also to have the same B. oleae population level. A diverse expression profile between the tolerant and susceptible genotype was clearly observed only in the presence of larvae. Within a species, the scale and quality of response to herbivores as well as the accumulation of defence compounds may vary significantly [47][48][49][50]. In Brassica oleracea, there was a very little overlap between transcriptional responses of two varieties with contrasting level of resistance to the cabbage aphid (Brevicoryne brassicae), underlying that intraspecific variation in susceptibility to insect pests can be also explained by differences in induced transcriptional changes [51]. The here described strong difference between the tolerant and the susceptible variety also implies a significant genotypespecific response in olive. This would be consistent with the frequently reported genetic diversity and the large diversified phenotypic traits of the plant and of olive oil products [52]. The data suggest a strong relevance of the genetic component in the determination of the tolerant phenotype [53,54] also because trees grew in the same environment using the same agronomic practices.
Plant defence pathways can be manipulated by also pathogens or pests [55,56]. For instance, glucose oxidase, one of the principal components of Helicoverpa zea saliva, suppresses some induced resistance in Nicotiana tabacum by directly inhibiting the woundsignalling molecule jasmonic acid and/or by antagonizing its interaction with other signalling pathways [57]. In tomato, the spider mite Tetranychus evansi is able to suppress the induction of genes involved in the induced plant defences, such as proteinase inhibitors [58]. Evidence for the defence suppression by herbivores is not as large as for defence induction yet, it is conceivable that plant resistance mechanisms promote the selection of counter-adaptive mechanisms in biotic stressors, especially in compatible interactions involving monophagous pests [59]. Our data suggest that the successful colonization of drupes by B. oleae is likely to be due to a weak reaction in the more susceptible variety and that drupes of the more tolerant variety are a less favourable nourishment because of a more active and composite molecular response [60].
The gene expression in the attacked drupes of cv. 'Ruveia' identified a number of sequences directly and indirectly involved in the induced resistance mechanism. Plant response to herbivorous pest encompasses a number of mechanisms [61] and therefore, several molecular processes are activated by the drupe defence to the fruit fly. Upregulated genes include those involved in the defence mechanisms against biotic stress, both wounding and pathogen attack, as well as abiotic stress, such as drought and high or low temperature. The overexpression of genes involved in different defense pathways suggests that the inducible response to B. oleae larvae has a degree of overlap with the response to pathogens, which was also noted for the interaction between olive and S. oleaginea [62]. Moreover, our data is consistent with other research works that indicated that beta-glucosidases are important element of the olive defense to B. oleae [14,15,37]. Beta-glucosidases promote the formation of toxic glutaraldehyde-like structure from oleuropein, a phenolic compound that has been linked to the different susceptibility of olive cultivars [16]. DEGs were also involved in the production or response to phytohormones and molecules involved in pest response (e.g., jasmonic acid and ROS).
Besides genes associated to stress resistance, the Gene Ontology analysis suggested that biological processes related to secondary metabolism, cation transport and transmembrane transport are also affected, giving reasons to believe that numerous changes in plant primary metabolism should occur in response to larva feeding [63]. Recently, it has been reported that B. oleae infestation causes significant changes in mineral elements (such as P, K, Fe and Mg) in fruits [23]. Also abiotic stress, such as a moderate drought, affects different metabolites during fruit development, including terpenes [64]. Collectively, the existing molecular and metabolomics data suggest that the defensive reaction of the tolerant cultivar to the persisting feeding of the B. oleae comprises a larger than anticipated metabolic reprogramming in infested tissues, yet to be fully described [7,65].

Conclusions
The B. oleae feeding influences pathways with a known role in defence, oxidative stress responses, and genes involved in plant structure and metabolism. Defence against the longlasting fruit fly larva feeding is a complex trait and involves multiple molecular mechanisms. The complexity of the drupe response suggests that a number of features, metabolites and signalling pathways effectively limit the fruit fly infestation. An additional level of complexity derives from the marked genotype-specific difference present in olive that is suggested by our study [66]. The extended gene expression difference between the tolerant and susceptible olive varieties under investigation also indicated that it would be possible to identify genetic factors that are associated with a higher level of tolerance to B. oleae.